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TE  SCATTERING  FROM  A  DIELECTRIC  COATED  CONDUCTING  STRIP: 


PROGRAM  “PBFSTRIP” 

I.  PBFSTRIP  Program  Description 

The  physical  basis  function  moment  method  implementation  was  coded  in  standard  FOR¬ 
TRAN.  The  program  was  written  to  make  maximal  use  of  COMMON  blocks  to  transfer  data  be¬ 
tween  program  units.  Capt  W.  Irvin  contributed  the  matrix  inversion  routine,  INVERT.  Frequency 
has  been  scaled  out  of  the  program,  so  all  length  variables  are  in  units  of  free-space  wavelength  A() 
and  lc0  =  2ir. 

The  program  performs  the  following  functions  sequentially. 

1 .  Define  elementary  constants,  read  geometry  and  material  parameters,  calculate  the  PBF  phase 
constants,  partition  the  dielectric  slab  into  rectangular  cells,  find  the  PBF  field  amplitudes, 
evaluate  the  surface  integrations  and  all  integrations  around  the  singularity. 

2.  Read  user-specified  match  points  and  fill  the  impedance  and  voltage  matrices. 

3.  Normalize  the  over-determined  system  of  equations. 

4.  Solve  the  resultant  3x3  system  of  equations  using  Gaussian  elimination.  Write  the  PB1 
amplitudes  to  an  output  file. 

5.  Calculate  residual  errors  from  the  least-squares  fit. 

6.  Calculate  and  write  equivalent  currents  along  the  strip  and  slab. 

7.  Calculate  the  bistatic  scattering  width. 

A  complete  listing  of  the  source  code  follows. 
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II.  Program  PBFSTRIP 


c********* ********** *************************************************** 


c  ** 

C  »»  PROGRAM  PBFSTRIP  ««  ** 
C  ** 
C  THIS  PROGRAM  CALCULATES  THE  SCATTERIIG  FROM  A  DIELECTRIC-  ** 
C  COATED  COIDUCTIIG  STRIP  USIIG  A  MOMEVT  METHOD  IMPLEMEKTATIOK  ** 
C  WITH  PHYSICAL  BASIS  FUICTIOIS  AID  LEAST-SqUARES  POUT  ** 
C  MATCHIHG.  ** 
C  ** 


C* ************************************************************* ******** 


c 

PROGRAM  PBFSTRIP 
C 

IMPLICIT  HONE 
INTEGER  N 

CHARACTER  DATESTRING*9 
C 

C** **************** **************************************************** 

C  SUBROUTINES  CALLED:  DEFINE,  DATE(INTRINSIC  FORTRAN), 

C  IMPEDANCE.MATRIX,  ERRORS,  FIHD.CURREHTS ,  RADIATE 

C  FUNCTIONS  CALLED:  NONE 

C  COMMON  BLOCKS:  ALL,  HAHKEL,  GAUSSIAN.QUADRATURE, 

C  FIELD.AMPLITUDES,  SURFACE. INTEGRALS,  IMPEDANCE, 

C  MATCH.POINT 

C 

C  »»  INTERNAL  VARIABLES  «« 

C  I  =  IMPLIED  DO-LOOP  INDEX 

C  DATESTRING  =  STRING  CONTAINING  CURRENT  DATE 

C 

C* ************************************************************ ******** 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DEFINE  THE  VARIABLES  IN  THE  COMMON  BLOCKS 

COMMON  BLOCK:  ALL  —  ELEMENTARY  VALUES 
PI  =  3.141592 . 

CJ  =  IMAGINARY  NUMBER  OF  UNIT  AMIPLITUDE 

EO  =  PERMITTIVITY  OF  FREE  SPACE 

MUO  =  PERMEABILITY  OF  FREE  SPACE 

ETA  =  IMPEDANCE  OF  FREE  SPACE 

KO  =  PHASE  CONSTANT  OF  FREE  SPACE,  EQUALS  2*PI 

THETA  =  ANGLE  OF  INCIDENCE,  IN  RADIANS 

ER  =  RELATIVE  PERMITTIVITY  OF  DIELECTRIC 

H  =  THICKNESS  OF  DIELECTRIC  SLAB 

V  =  WIDTH  OF  SLAB  AND  STRIP 

F  =3X1  VECTOR  CONTAINING  Y  PHASE  CONSTANTS 

G  =  3X1  VECTOR  CONTAINING  X  PHASE  CONSTANTS 

XNODES=  NUMBER  OF  NODES  IN  THE  X  DIRECTION 
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C  Y10DES=  NUMBER  OF  IODES  II  TEE  Y  DIRECTIOI 
C  DELX  =  WIDTH  OF  IITEGRATIOI  CELL 

C  DELY  =  HEIGHT  OF  IHTEGRATIOI  CELL 

C 

C  COMMON  BLOCK:  HANKEL  —  VALUES  FOR  HANKEL  FUNCTIONS  APPROXIMATIONS 
C  JO,  YO,  MAGO,  PHASEO,  ALPHA,  BETA  = 

C  CONSTANTS  USED  IN  FUNCTION  HANKO 

C  Jl,  Yl,  MAGI,  PHASE1  =  CONSTANTS  USED  II  FUNCTION  HANK1 

C 

C  COMMON  BLOCK:  GAUSSIAI.QUADRATURE  —  IODES  AID  WEIGHTS 
C  N8  =  8-POINT  GAUSSIAN  QUADRATURE  IODES 

C  W8  =  8-POINT  GAUSSIAI  QUADRATURE  WEIGHTS 

C  14  =  4-POINT  GAUSSIAI  QUADRATURE  IODES 

C  W4  =  4-POINT  GAUSSIAI  QUADRATURE  WEIGHTS 

C 

C  COMMON  BLOCK:  FIELD.AMPLITUDES  —  PBF  FIELD  AMPLITUDES 
C  CX  =  3X1  VECTOR  OF  MULTIPLIERS  OF  X  COMPOIEITS  OF  PBF* 8 

C  CY  =  3X1  VECTOR  OF  MULTIPLIERS  OF  Y  COMPONENTS  OF  PBF’s 

C  CZ  =  3X1  VECTOR  OF  MULTIPLIERS  OF  Z  COMPOIEITS  OF  PBF’s 

C 

C  COMMON  BLOCK:  SURFACE_IITEGRALS  —  PRE-CALCULATED  SURFACE  INTEGRALS 
C  AID  INTEGRALS  ABOUT  THE  SINGULARITY 

C  SURFINT(I, J,I,K)  =  MATRIX  CONTAIIIIG  SURFACE  IITEGRATIOI  OVER 
C  DIELECTRIC  CELL  WITH  LOWER  RIGHT  CORNER  AT  (I*DELX , J*DELY) , 

C  WHERE  INTEGRAND  INVOLVES  Ith  PBF  AID  EITHER  COSIIE  (K=l)  OR 

C  SINE  (K=2) . 

C  SIIGINT(K.I)  =  MATRIX  CONTAIIIIG  IITEGRATIOIS  ABOUT  THE 
C  SINGULARITY.  SEE  SUBROUTINE  SINGULAR.IITEGRAL  FOR  MEANING 

C  OF  K  AND  I. 

C 

C  COHMON  BLOCK:  IMPEDANCE  —  SYSTEM  OF  EQUATIONS 
C  Z(M,I)  =  MATRIX  CONTAINING  UINORMALIZED  IMPEDANCE  MATRIX 

C  ELEMENT  FOR  Mth  MATCH  POINT  AID  Ith  PBF. 

C  V(M)  =  VECTOR  CONTAINING  UIIORMALIZED  VOLTAGE  MATRIX  ELEMENT 

C  FOR  Mth  MATCH  POINT. 

C  ZN(M,N)  =  3X3  MATRIX  COITAIIIIG  NORMALIZED  IMPEDANCE  MATRIX 

C  VI (N)  =  VECTOR  COITAIIIIG  NORMALIZED  VOLTAGE  MATRIX 

C  CUR(N)  =  VECTOR  COITAIIIIG  PBF  AMPLITUDES 

C  CN  CONDITION  NUMBER  OF  ZI 

C  IUMMPS  =  TOTAL  I  UMBER  OF  MATCH  POUTS 

C 

C  COMMON  BLOCK:  MATCH.POIIT  —  VALUES  ASSOCIATED  WITH  CURRENT  MATCH 
C  POUT 

C  I  =  MATCH  POUT  X-COORDIIATE  IIDEX 

C  J  MATCH  POUT  Y-COORDIIATE  IIDEX 

C  X  =  MATCH  POUT  X-COORDIIATE  (  =  I*DELX) 

C  Y  MATCH  POUT  Y-COORDIIATE  (  =  J*DELY) 

C  H02PY  =  H  OVER  2,  PLUS  Y 
C  H02PY2=  H02PY  SQUARED 

C  H02MY  =  B  OVER  2,  MIIUS  Y 
C  H02MY2=  H02MY  SQUARED 


3 


C  V02PX  =  V  OVER  2,  PLUS  X 

C  V02PX2*  W02PX  SQUARED 

C  V02NX  =  V  OVER  2.  MIVUS  X 

C  V02MX2=  V02MX  SQUARED 

C 

C - CONNOR  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  NUO, ETA, KO, THETA, ER,H,V,F(3) ,G(3) ,DELX,DELY 
INTEGER  XIODES,  YIODES 

CONNOR  /  ALL  /  PI,  CJ,  EO,  NUO,  ETA,  KO,  THETA,  ER,  H,  W,  F,  G, 

*  XIODES ,  YIODES,  DELX,  DELY 
C 

REAL*8  J0(0:6) ,  Y0(0:6),  BETA(0:6),  NAG0(0:6),  PHASE0(0:6) 
COMPLEX* 16  ALPHA(0:6) 

REAL*8  Jl(0:6) ,  Yl(0:6),  NAG1(0:6),  PHASE1(0:6) 

CONNOR  /  HARKEL  /  JO,  YO,  ALPHA,  BETA,  NAGO,  PHASEO, 

*  Jl,  Yl,  NAG1 ,  PHASE1 
C 

REAL* 8  V8(l:8),  18(1:8),  W4(l:4),  14(1:4) 

CONNOR  /  GAUSSIAR.qUADRATURE  /  W8,  R8,  w4,  n4 
C 

COMPLEX* 16  CX(3)  ,  CY(3),  CZ(3) 

CONNOR  /  FIELD.ANPLITUDES  /  CX,  CY,  CZ 
C 

C0NPLEX*16  SURFIRT(-49:50,-S:S, 1 :3, 1 :2) ,  SIRGIRT(4,3) 

CONNOR  /  SURFACE. INTEGRALS  /  SURFIRT,  SIRGIRT 
C 

INTEGER  RUMMPS 
REAL* 8  CR 

COMPLEX* 16  Z(40,3) ,  V(40) ,  ZR(3,3),  VR(3),  CUR(3) 

CONNOR  /  IMPEDANCE  /  Z,  V  ,  ZR,  VR,  CUR,  CR ,  RUMMPS 
C 

INTEGER  I,  J 

REAL*8  X,  Y,  H02PY,  H02PY2,  H02NY,  H02MY2, 
ft  V02PX ,  V02PX2,  V02NX ,  V02MX2 

CONNOR  /  MATCH.POIRT  /  I,  J,  X,  Y,  H02PY,  H02PY2,  H02NY ,  H02MY2, 
ft  V02PX,  V02PX2 ,  V02NX,  V02NX2 

C 

C  OPEN  FILES  FOR  INPUT  AND  OUTPUT 
C 

OPEN  (URIT=19 ,  FILE* ’ PBFSTRIP . II  ’,  STATUS* ’ OLD ’ ) 

OPEN  (URIT=20,  FILE* ’PBFSTRIP. OUT’,  STATUS* 'NEW* ) 

OPEN  (URIT=31,  FILE=’PBFSTRIP_CUR.PL1’ ,  STATUS* ' REV ’ , 
ft  FORM* 'UNFORMATTED’) 

OPEN  (URIT=22 ,  FILE* ’ PBFSTRIP. FLD’ ,  STATUS* ’ NEW ’ ) 

OPEN  (UIIT=32,  FILE*’PBFSTRIP_FLD.PL1 ' ,  STATUS* ’ NEW ’ , 
ft  FORM* ’ UNFORMATTED ’ ) 

CALL  DEFIIE 
C 

C  WRITE  HEADER  INFORMATION  TO  THE  OUTPUT  FILES 
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c 

CALL  DATE  (DATESTRIIG) 

WRITE  (20,98)  ER,  DATESTRIIG,  W,  H,  MODES,  MODES, 
ft  THETA*180/PI , (G(M) ,  F(M) ,  1=1,3) 

98  FORMAT  ('  Relative  Permittivity  =  »,  F10.4,  10X,  A9/ 

ft  »  Slab  Width  =  ',  F10.4, 

ft  ‘  Slab  Thickness  =  »,  F10.4/ 

ft  *  Modes  (horizontal)  =  ’,  13,  *  (vertical)  =’,  12/ 

ft  ’  Angle  from  lormal  =  ’ ,  F10.4,  '  deg'/ 

ft  '  Phase  Constants:  ',  3(' (' ,F7.3, ' , » ,F7.3, ')  ')  ) 

WRITE  (21,99)  ER,  DATESTRIIG,  W,  H,  THETA* 180/PI 
WRITE  (22,99)  ER,  DATESTRIIG,  W,  H,  THETA*180/PI 

99  FORMAT  (’  Relative  Permittivity  =  ’ ,  F10.4,  10X,  A9/ 

ft  *  Slab  Width  =  F7.2,  *  Slab  Thickness  =  ',  F7.3, 

ft  ’  Angle  from  lormal  =  ' ,  F7 . 3  ) 

C 

C  CALCULATE  THE  PHYSICAL  BASIS  FUICTIOI  AMPLITUDES 
C 

CALL  IMPED AICE_MATRIX 
C 

C  CALCULATE  THE  MEAI-SQUARE  ERRORS,  EQUIVALEIT  CURREITS,  AID  FAR-ZONE 
C  SCATTERED  FIELDS. 

C 

CALL  ERRORS 
CALL  FIID_CURREITS 
CALL  RADIATE 
C 

C  CLOSE  IIPUT  AMD  OUTPUT  FILES 
C 

CLOSE  (19) 

CLOSE  (21) 

CLOSE  (22) 

CLOSE  (31) 

CLOSE  (32) 

C 

STOP 

EID 
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2.1  Subroutine  DEFINE 


C****«***************************************************************** 


C  ** 
C  »»  SUBROUTIIE  DEFIIE  ««  ** 
C  ** 
C  THIS  SUBROUTIIE  IIITIALI2ES  THE  COISTAITS  II  COKMOH  ** 
C  BLOCKS  ALL,  HAIKEL,  GAUSSI AB .QUADRATURE ,  *♦ 
C  FIELD.AMPLITUDES ,  AID  SURFACE.IITEGRALS .  ** 
C  ** 


C********************************************************************** 

C  CALLED  BY:  MAIN 

C  SUBROUTIIES  CALLED:  PHASE.COISTAITS ,  MATRIX .FILL, 

C  DEFIIE.IMIABC 

C  FUICTIOIS  CALLED:  IOIE 

C  COMMOH  BLOCKS:  ALL,  HAIKEL,  GAUSSIAI.QUADRATURE , 

C  FIELD. AMPLITUDES,  SURFACE.IITEGRALS 

C 

C  »»  UTERI AL  VARIABLES  «« 

C  THETA.DEG  =  AIGLE  OF  IICIDEICE,  II  DEGREES 
C  SAMPLE  =  I UMBER  OF  BODES  PER  FREE-SPACE  WAVELENGTH 
C  I  =  DO-LOOP  IIDEX 

C 

C  »»  DATA  IIPUT  FROM  CALLIIG  ROUTIIE  «« 

C  IOIE 

C 

C  »»  DATA  OUTPUT  «« 

C  ALL  DATA  IS  OUTPUT  VIA  THE  COMMOI  BLOCKS 
C 

C********************************************************************** 
SUBROUTIIE  DEFIIE 
C 

IMPLICIT  IOIE 
REAL* 8  THETA.DEG,  SAMPLE 
IITEGER  I 
C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL* 8  PI,  EO,  MUO, ETA, KO, THETA. ER,H,W,F(3) ,G(3) ,DELX ,DELY 
IITEGER  XIODES,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUC,  ETA,  KO,  THETA,  ER,  H,  W.  F,  G, 

*  XIODES,  YIODES,  DELX,  DELY 

C 

REAL* 8  J0(0:6) ,  Y0(0:8),  BETA(0:6),  MAG0(0:6),  PHASE0(0:6) 
COMPLEX* 16  ALPHA (0:6) 

REAL*8  Jl(0:6) ,  Yl(0:8),  MAG1(0:6),  PHASE1(0:6) 

COMMOI  /  HAIKEL  /  JO,  YO,  ALPHA,  BETA,  MAGO,  PHASEO, 
t  Jl.  Yl,  MAGI,  PHASE1 

C 
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REAL* 8  W8(l :8) ,  18(1:8),  W4(l:4),  14(1:4) 

CONNOR  /  GAUSS I AN. QUADRATURE  /  W8,  18,  s4,  n4 
C 

CONPLEX*16  CX(3) ,  CY(3) ,  CZ(3) 

CONNOR  /  FIELD. ANPLITUDES  /  Cl,  CY,  CZ 
C 

C0NPLEX*16  SURFIRT (-49 :  50 ,  -5 : 5 , 1 : 3 , 1 : 21 ,  SIRGIRT(4,3) 
CONNOR  /  SURFACE. IRTEGRALS  /  SURFIRT,  SIRG1RT 
C 
C 

~  DEFIRE  SONE  ELENERTARY  CORSTARTS 
C 

PI=3. 14159  26S35  89793  DO 
CJ=(O.DO,  1 .DO) 

E0=  8.85418  53368D-12 
K0=2.D0*PI 

NU0=1. 25663  70614D-06 
ETA=DSQRT(NUO/EO) 

C 

C  READ  IR  CASE-SPECIFIC  PARAMETERS 
C 

READ  (19,56)  THETA.DEG ,  ER,  H.  W,  SANPLE 
56  FORNAT  (T30,  F15.8) 

C 

C  CONVERT  THETA.DEG  TO  RADIANS 
C 

THETA=THETA_DEG*PI/ 180 . DO 
C 

C  DEFINE  NODES  SUCH  THAT  THE  SEPARATION  BETWEEN  THEN  IS 
C  LESS  THAN  ORE  WAVELENGTH  IH  THE  DIELECTRIC.  ENSURE  AT  LEAST 
C  ONE  LAYER  OF  NODES  OH  THE  CONDUCTOR  AND  ONE  IN  THE  DIELECTRIC. 

C 

XNODES=IHT(NP_PER_WAVELERGTH*W) - 1 
YHODES=INT (NP_PER_WAVELENGTH*H) 

IF  (YNODES.LT. 2)  YN0DES=2 
C 

C  CALCULATE  DELX  AND  DELY  BASED  OR  XRODES  AND  YRODES .  ENSURE  THAT 
C  DELY  IS  LESS  THAN  DELX. 

C 

DELX=W/ (XRODES+1) 

10  DELY=H/YNODES 

IF  (DELY. GE. DELX)  THEN 
YNODES=YNODES+ 1 
GOTO  10 
ERDIF 
C 

C  CHECK  IF  TOO  NARY  NODES  FOR  THE  DIHENSIONALITY  OF  SURFINT. 

C 

IF  (XRODES. GT. 60)  THEN 

WRITE  (6,*)  'Strip  too  wide  in  subroutine  DEFINE' 

STOP 
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ELSEIF  (YIODES.GT.S)  TFEI 

WRITE  (6,*)  'Dielectric  slab  too  thick  in  subroutine  DEFINE’ 
STOP 
EMDIF 
C 

C  DEFINE  THE  PHASE  CONSTANTS  F  AND  G 
C 

CALL  PHASE.CONSTANTS 
C 

C  JO,  YO,  MAGO,  PHASE0.J1,  Yl,  MAGI,  AND  PHASE1  ARE  DEFINED  IN  THE 
C  HANKEL_COEFFICIENTS  BLOCK  DATA  SUBROUTINE.  CALCULATE  ALPHA  AND 
C  BETA  BASED  ON  THEM. 

C 

DO  15  1=0,6 

BETA(I)=2.D0*J0(I)/PI 

15  ALPHA(I)=J0(I)+CJ*(BETA(I)*DL0G(2 .DO)-YO(I)) 

C 

C  GAUSSIAN  QUADRATURE  WEIGHTS  AND  NODES  DEFINED  IN 
C  GAUSS I AN_QUADRATURE  BLOCK  DATA  SUBROUTINE. 

C  DEFINE  THE  FIELD  AMPLITUDES,  CX,  CY,  AND  CZ 
C 

CX ( 1 ) =C J*DSQRT ( 1 . DO-DSIN (THETA) **2/ER) 

CY ( 1 ) =DSIN (THETA) /DSQRT(ER) 

CZ ( 1 ) =DSQRT(ER*EO/MUO) 

CX (2) =F( 2) ♦ETA/ (CJ*ER*KO) 

CY(2)=-G(2)*ETA/(K0*ER) 

CZ(2)=-1 .DO 

CX(3)=F(3)*ETA/ (CJ*ER*KO) 

CY(3)=-G(3)*ETA/(K0*ER) 

CZ(3)=-1 .DO 
C 

C  CALCULATE  THE  ELEMENTS  OF  SURFINT. 

C 

CALL  MATRIX.FILL 
C 

C  CALL  SINGULAR.INTEGRAL  TO  CALCULATE  THE  INTEGRATIONS  AROUND  THE 
C  SINGULAR  POINT. 

C 

CALL  SINGULAR_INTEGRAL 
C 

RETURN 

END 
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2.1.1  Block  Data  Subroutine  HANKEL.COEFFICIENTS 


C**** ****************************************************************** 


c  ** 

C  »»  BLOCK  DATA  SUBROUTINE  HANKEL.COEFFICIENTS  ««  ** 
C  ** 
C  THIS  SUBROUTIHE  INITIALIZES  COISTAITS  USED  II  FUICTIOI  ** 
C  HAHKO.  VALUES  COME  FROM  AMS-5S,  EQUATIONS  9.4.1  THROUGH  ** 
C  9.4.6.  ** 
C  ** 


C******* *************************************************************** 

BLOCK  DATA  HANKEL.COEFFICIENTS 

REAL*8  J0(0:6) ,  Y0(0:6),  BETA(0:6),  MAG0(0:6),  PHASE0(0:6) 
COMPLEX* 16  ALPHA (0:6) 

REAL* 8  Jl(0:6),  Yl(0:6),  MAG1(0:6),  PHASE1(0:6) 

COMMON  /  HANKEL  /  JO,  YO,  ALPHA,  BETA,  MAGO,  PHASEO, 
k  Jl,  Yl,  MAGI,  PHASE1 

DATA  (JO(I),  1=0,6)  /  + 1 . OOOOOOODO , -2 . 2499997D0 , +1 . 2656208D0 , 
k  -0 . 3163866D0 , +0 . 0444479DO ,-0 . 0039444D0 , +0 . 0002100DO/ 

DATA  (YO(I) ,  1=0,6)  /+0.36746691D0,+0.60559366D0,-0.74350384D0, 
k  +0 . 25300 1 1 7D0 , -0 . 0426 1 214D0 , +0 . 004279 16D0 , -0 . 00024846D0/ 

DATA  (MAGO (I) ,  1=0,6)  /+O.79788456D0,-0.00000077D0,-0.O0552740D0, 
ft  -0 . 000095 12D0 ,+0 .00 137237D0 ,-0 .00072805DO, +0 . 00014476D0/ 

DATA(PHASEO(I) ,  1=0 . 6)/-0 . 78539816D0 ,-0 . 04166397D0 , -0 . 00003954D0 , 
ft  +0 . 00262573D0 . -0 . 00054125D0 , -0 . 00029333D0 , +0 . 00013558D0/ 

DATA  (J1(I),  1=0,6)  /  +0 . 5000000CD0 , -0 . 56249985D0 , +0 . 21093573D0 , 
ft  -0.03954289D0.+0. 00443319D0, -0.0003 1761D0, +0. 00001 109D0/ 

DATA  (Y1(I)  ,  1=0,6)  /  -0 . 6366198D0 , +0 . 221 209 IDO , +2 . 1682709D0 , 
ft  -1 . 3164827D0 , +0 . 3123951D0 ,-0 . 0400976DO , +0 . 0027873D0/ 

DATA  (MAGl(I) ,  1=0,6)/  +0.79788456DO,+0.00000156DO,+0.01659667DO, 
ft  +0 . 000171 05D0 , -0 . 0024961 IDO , +0 . 001 13653D0 , -0 . 00020033D0/ 

DATA(PHASE1(I),  I=0,6)/-2.35619449DO,+0.12499612DO,+0.000''5650DO, 
ft  -0 . 00637879D0 , +0 . 00074348D0 ,+0 . 00079824D0 , -0 . 00029166D0/ 

END 
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2.1.2  Block  Data  Subroutine  GA  USS1A  N.Q  UA  DRA  TURE 


C* ********************************************************************* 


c  ** 

C  »»  BLOCK  DATA  SUBROUTIIE  GAUSSIAI_QUADRATURE  ««  ** 
C  ** 
C  THIS  SUBROUTIIE  IIITIALIZES  WEIGHTS  AID  IODF.S  FOR  8-POIIT  ** 
C  GAUSSIAI  QUADRATURE  IITEGRATIOI .  VALUES  CONE  FROM  AMS-55,  ** 
C  TABLE  25.4  ** 
C  ** 


C*** ****** ************************* ******* ************ ***************** 
BLOCK  DATA  GAUSSIAI.QUADRATURE 
REAL*8  V8(l:8) ,  18(1:8),  W4(l:4),  14(1:4) 

COHMOI  /  GAUSSIAI.QUADRATURE  /  W8,  18,  W4,  14 
DATA  (V8(I),  1=1,8)  / 

A  0.10122  85362  90376D0,  0.22238  10344  53374D0, 

A  0.31370  66458  77887D0,  0.36268  37833  78362D0, 

A  0.36268  37833  78362D0,  0.31370  66458  77887D0, 

A  0.22238  10344  53374D0,  0.10122  85362  90376D0/ 

DATA  (18(1),  1=1,8)  / 

A  -0.96028  98564  97536D0, -0.79666  64774  13627D0, 

A  -0.62553  24099  16329D0.-0. 18343  46424  95?  '0, 

A  0.18343  46424  95650D0,  0.52553  24099  161  .„0, 

A  0.79666  64774  13627D0,  0.96028  98564  97536D0/ 

DATA  (W4(I) ,  1=1,4)  / 

A  0.34785  48451  37454D0,  0.65214  51648  62546D0, 

A  0.65214  61648  62646D0,  0.34785  48451  37454D0/ 

DATA  (14(1),  1=1,4)  / 

A  -0.86113  63115  94053D0, -0.33998  10435  84856D0, 

A  0.33998  10435  84856D0,  0.86113  63115  94053D0/ 

END 
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2.1.S  Subroutine  PHASE.CONSTANTS 


C* ************************************************************** ******* 

c  ** 

C  »»  SUBROUTINE  PHASE.CONSTANTS  ««  ** 

C  ** 

C  THIS  SUBROUTINE  CALCULATES  THE  PHASE  COISTA1TS  OF  THE  ** 

C  THREE  PHYSICAL  BASIS  FU1CTIOIS.  F  IS  THE  PHASE  CONSTANT  ** 

C  II  THE  Y  DIRECTIOI  AID  G  IS  THE  PHASE  COISTAIT  II  THE  ** 

C  X  DIRECTIOI.  1=1  IS  THE  FORCED  WAVE,  1=2  IS  THE  FORWARD  ** 

C  SURFACE  WAVE,  AID  1=3  IS  THE  REVERSE  SURFACE  WAVE.  THE  ** 

C  IEWTOI-RAPHSOI  ROOT-FIIDIIG  ALGORITHM  IS  USED  TO  FIID  ** 

C  F(2) ,  WHICH  IS  THE  ROOT  OF  THE  TRAISCEIDEITAL  EqUATIOI  ♦* 

C  F(2)*TAN(F(2)*H)-ER*SqRT((ER-l)*K0**2-F(2)**2)=0,  WHERE  ** 

C  F(2)  LIES  II  THE  IITERVAL  (0,  PI/(2*H)).  ** 

C  ** 

c********************************************************************** 
C  CALLED  BY:  DEFIIE 

C  SUBROUTIIES  CALLED:  IOIE 

C  FUICTIOIS  CALLED:  IOIE 

C  COMMOI  BLOCKS:  ALL 

C 

C  »»  IITERIAL  VARIABLES  «« 

C  TOL  =  THRESHOLD  FOR  COIVERGEICE  OF  IEWTOI-RAPHSOI  ALGORITHM 

C  SLOPE  =  DERIVATIVE  OF  TRAISCEIDEITAL  EqUATIOI 

C  OLDX  =  OLD  GUESS  FOR  F(2),  USED  TO  FIID  IEW  GUESS  FOR  F(2) 

C  OLDY  =  VALUE  OF  TRAISCEIDEITAL  EqUATIOI  EVALUATED  AT  OLDX 

C  I  ITERATION  COUNTER  IN  IEWTON-RAPHSON  ALGORITHM 

C 

C  »»  DATA  OUTPUT  «« 

C  F  AND  G  VECTORS  RETURNED  VIA  COMMOI  BLOCK  "ALL" 

C*** ******************************************************* ************ 
SUBROUTINE  PHASE.COISTAITS 
C 

IMPLICIT  IOHE 

REAL* 8  TOL,  OLDX,  OLDY,  SLOPE 
INTEGER  I 
C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MUO, ETA, KO, THETA, ER,H,W,F(3),G(3),DELX,DELY 
IITEGER  XIODES,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  W,  F,  G, 
ft  XIODES,  YIODES,  DELX,  DELY 

C 
C 

C  DEFIIE  THE  PHASE  CONSTANTS  OF  THE  FORCED  WAVE 
C 

F(l)  =  K0*DSqRT(ER-DSII(THETA)**2) 

G(l)  =  KO*DSII(THETA) 
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c 

C  DEFINE  TOL,  WHICH  CONTROLS  NEWTON-RAPHSON  C01VERGE1CE.  IF  ER  IS 
C  VERY  SMALL,  WE  CAI  DEFIIE  F  AID  G  EXPLICITLY. 

C 

TOL  =  l.E-08 
IF  (ER-l.LT.TOL)  THEI 
F(2)=0.0 
F(3)=0.0 
G(2)=K0 
G(3)=-K0 
RETURI 
EIDIF 
C 

C  MAKE  A  FIRST  GUESS  AT  F(2)  BASED  01  ITS  MAXIMUM  POSSIBLE  VALUE 
C 

F(2)=MII(DSQRT(ER-1 . )*K0 ,  0 . 99*PI/(2 . *H)  ) 

C 

C  USE  THE  IEWTOH-RAPHSOH  METHOD  TO  ITERATIVELY  FIID  F(2) . 

C  INITIALIZE  COUNTER  I. 

C 

1=0 

10  1=1+1 
C 

C  TRANSITION  THE  CURRENT  GUESS  TO  OLDX,  AND  EVALUATE  THE 
C  TRANSCENDENTAL  EQUATION  AS  OLDY.  ALSO,  FIND  THE  DERIVATIVE 
C  OF  THE  TRANSCENDENTAL  EQUATION  AS  SLOPE. 

C 

OLDX  =  F(2) 

OLDY  =  OLDX*OLDX+TAN(OLDX*H)+TAN(OLDX*H)  - 
*  ER*ER*( (ER-1 . )*KO*KO-OLDX*OLDX) 

SLOPE  =  2*0LDX*TAN(0LDX*H) *TAN (OLDX*H) 
t  +  2*0LDX*0LDX*H*TAN(0LDX*H)/C0S(0LDX*H)**2 
A  +  2*ER*ER*0LDX 

C 

C  CALCULATE  THE  NEW  GUESS,  F(2) ,  ACCORDING  TO  NEWTOI-RAPHSON 
C 

F(2)  =  MIN (OLDX  -  OLDY/SLOPE,  K0*DSQRT(ER-1 . )) 

C 

C  IF  THE  NEW  GUESS  LIES  OUTSIDE  THE  POSSIBLE  RANGE  OF  VALUES,  THEN 
C  SOMETHING  IS  DREADFULLY  WRONG.  RETURN  AID  FLAG  THE  ERROR. 

C 

IF  ( (F(2) .LE.O. ) .OR. (F(2) .GT.PI/(2 . *H) ) )  THEN 
WRITE  (6,*)  'ERROR  IN  SUBROUTINE  "COEF"’ 

RETURN 

ENDIF 

C 

C  TEST  FOR  CONVERGENCE.  IF  CONVERGENCE  IS  NOT  REACHED,  ITERATE  AGAIN. 
C  ONLY  ITERATE  60  TIMES 
C 

IF  ( (DABS ((F( 2) -OLDX) /OLDX) .GT.TOL) . AID . (I .LT. 60) )  GOTO  10 
C 
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C  CHECK  TO  SEE  IF  MAXIMUM  IUMBER  OF  ITERATIOIS  HAS  BEEI  EXCEEDED 
C 

IF  (I.GE.50)  THEI 

WRITE  (6,*)  'Ion-convergence  in  SUBROUTINE  COEF' 

RETURN 

EIDIF 

C 

C  GENERATE  F(3),  G(2),  AID  G(3)  BASED  ON  F(2) 

C 

F(3)  =  F(2) 

G(2)  =  DSQRT(ER*KO*KO  -  F(2)*F(2)) 

G(3)  =  -G(2) 

C 

RETURI 

END 
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2.1.4  Subroutine  MATRIX-FILL 


O********************************************************************* 


c  ** 

C  »»  SUBROUTINE  MATRIX.FILL  ««  ** 
C  ** 
C  THIS  SUBROUTIIE  FILLS  THE  SURFIIT  MATRIX.  EACH  ELEMEHT  ** 
C  COITAIVS  THE  IfTEGRATIOV  OF  A  FUICTIOV  OVER  A  RECTAVGULAR  ** 
C  CELL  OF  DIMEISIO*  DELX-by-DELY.  THE  LOWER-RIGHT  CORHER  OF  ** 
C  THE  CELL  IS  LOCATED  AT  (I*DELX,  J*DELY) .  THE  FUICTIOI  ** 
C  IITEGRATED  DEPENDS  OH  I  AID  K.  ** 
C  ** 


C* ********************************************************************* 

C  »»  CALLED  BY:  DEFIIE 

C  »»  SUBROUTIIES  CALLED:  MATRIX.ELEMENT,  CORKERS 

C  »»  FUHCTIOIS  CALLED:  IOIE 

C  »»  COMMOI  BLOCKS  USED:  ALL,  SURFACE.IITEGRALS 
C 

C  »»  DATA  FROM  MATRIX_ELEMENT  «« 

C  SURFIIT(I , J, 1 :3,1 : 2)  WHERE  (I.J)  DEFIIE  A  CELL  IOT  TOUCHIMG 
C  THE  SIIGULAR  POUT 

C 

C  »»  DATA  FROM  CORIERS  «« 

C  SURFINT(I, J,  1 :3, 1 :2)  WHERE  (I.J)  DEFIIE  A  CELL  TOUCHIKG 

C  THE  SIIGULAR  POUT 

C 

C  »»  IITERIAL  VARIABLES  «« 

C  I,  J  =  DEFIIE  COORDIIATES  OF  LOWER-RIGHT  CORIER  OF 

C  IITEGRATIOI  CELL 

C  I  PHYSICAL  BASIS  FUICTIOI  IIDEX 

C  K  IITEGRATIOI  FUICTIOI  IIDEX,  DEFIIES  TRIG.TERM 

C 

C  »»  DATA  OUTPUT  TO  CALLIIG  ROUTIIE  «« 

C  IIT  =  MATRIX  COITAIIIKG  IITEGRATIOIS  OF  (I.J)  CELL  FOR  THREE 
C  PHYSICAL  BASIS  FUICTIOIS  AID  TWO  TRIGONOMETRIC  FUICTIOIS 

C 

C********************************************************************** 

SUBROUTIIE  MATRIX.FILL 
C 

IMPLICIT  IOIE 
IITEGER  I,  J,  I 
C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MUO, ETA, KO, THETA, ER,B,W,F(3) ,G(3) .DELX.DELY 
IITEGER  XIODES,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA.  KO,  THETA,  ER,  H,  W,  F,  G, 
t  XIODES,  YIODES,  DELX,  DELY 

C 

COMPLEX* 16  SURFIIT(-49:60,-5:6,l:3,l:2),  SIIGIIT(4,3) 
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COMMOI  /  SURFACE. IITEGEALS  /  SURFIIT,  SIIGIKT 
C 

C  CALCULATE  THE  ELEMEHTS  WITH  CELLS  TOUCHIIG  THE  SIIGULAR  POINT 
C 

CALL  CORIERS 
C 

C  FILL  ELEMEHTS  WITH  CELLS  HOT  TOUCHIIG  THE  SIIGULAR  POUT. 

C  START  THE  I  AID  J  LOOPS 
C 

DO  10  I  =  1-XIODES,  XIODES 
DO  10  J  =  0,  YHODES-1 

CALL  MATRIX_ELEMEHT(I,  3 ) 

C 

C  EXPLOIT  SYMMETRY  BETWEEI  THE  (I,J,H,K)  AID  (I,-1-J,I,K)  ELEMEHTS 
C 

DO  10  H=1 ,3 

SURFIHT (I,-1-J,H,1) =-SURFIHT ( I , J , I , 1 ) 

SURFIHT(I ,-1-J ,H,2)=  SURFIIT(I, J,H,2) 

10  COHTIHUE 

C 

RETURI 

EHD 
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2.1.11  Subroutine  MATRIX-ELEMENT 


C********************************************************************** 


C  ** 
C  »»  SUBROUT  I  IE  KATRIX.ELEHENT  ««  ** 
C  ** 
C  THIS  SUBROUTIIE  CALCULATES  THE  ELEMEITS  OF  ** 
C  SURFIIT(I, J,*,K)  THAT  ARE  DEFIIED  BY  THE  CELL  REFERENCE  ** 
C  POIIT  COORDIIATES  (I,J).  THE  CELL  CAIIOT  DIRECTLY  TOUCH  THE  ** 
C  SI1GULAR  POIIT.  IITEGRATIOI  IS  DOIE  USIIG  8-POIIT  ** 
C  GAUSSIAI  QUADRATURE.  ** 
C  ** 


C** *************************************** ***************************** 

C  »»  CALLED  BY:  MATRIX.FILL 

C  »»  SUBROUTIIES  CALLED:  IOIE 

C  »»  FUICTIOIS  CALLED:  HAIKO 

C  »»  COMMON  BLOCKS:  ALL,  GAUSSIAN.QUADRATURE ,  SURFACE.INTEGRALS 
C 

C  »»  INTERNAL  VARIABLES  «« 

C  U,  V  =  COORDINATES  OF  POINT  AT  WHICH  FUNCTION  IS  EVALUATED 

C  DURING  GAUSSIAN  QUADRATURE 

C  TRIG.TERM  =  PART  OF  FUNCTION  INTEGRATED 

C  HANK.TERM  =  PART  OF  FUNCTION  INTEGRATED 

C  EXP_TERM  =  PART  OF  FUNCTION  INTEGRATED 

C  TERM  =  TEMPORARY  STORAGE  MATRIX  OF  INNER  INTEGRATIONS 

C  N  PHYSICAL  BASIS  FUNCTION  INDEX 

C  K  *  INTEGRATION  FUNCTION  INDEX,  DEFINES  TRIG.TERM 

C  II  INNER  GAUSSIAN  QUADRATURE  INDEX 

C  JJ  OUTER  GAUSSIAN  QUADRATURE  INDEX 

C 

C  »»  DATA  INPUT  FROM  CALLING  ROUTINE  «« 

C  I ,  J  =  COORDINATES  OF  RECTANGULAR  CELL 

C 

C  »»  DATA  OUTPUT  TO  CALLING  ROUTINE  «« 

C  SURFINT(I, J,N,K)  WHERE  (I,J)  DEFINE  A  CELL  NOT  TOUCHING 

C  THE  SINGULAR  POINT 
C 

C********************************************************************** 
SUBROUTINE  MATRIX.ELEMENT  (I.  J) 

C 

IMPLICIT  NONE 

REAL* 8  U,  V,  TRIG.TERM 

COMPLEX* 16  HANKO,  TERM(1 :3 , 1 : 2) ,  EXP.TERM,  HANK.TERM 
INTEGER  I.  J,  N,  K,  II,  JJ 
C 

C - COMMON  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  KUO, ETA, KO, THETA, ER,H,W,F(3) ,G(3) .DELX.DELY 
INTEGER  XNODES ,  YNODF.S 

COMMON  /  ALL  /  PI,  CJ,  EO,  KUO,  ETA,  KO,  THETA,  ER,  H,  W,  F,  G, 
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c 


XIODES,  YIODES,  DELX.  OELY 


* 


C 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


REAL*8  V8(l:8),  18(1:8),  V4(l:4),  14(1:4) 

CONNOI  /  GAUSSIAI.QUADRATURE  /  V8,  18,  s4,  n4 

C0MPLEX*16  SURFIIT(-49:S0,-5:5,1:3,1:2),  SIIGIIT(4,3) 
COMHOS  /  SURFACE.IITEGRALS  /  SURFIIT,  SIIGIIT 

IF  (I,J)  DEFIIE  A  CELL  TOUCHIIG  THE  SIIGULAR  POUT,  RETURI 

IF  (  ((I.EQ.O) .AID. (J.EQ.  0))  .OR. 

*  ((I.EQ.O). AID. (J.EQ.-l))  .OR. 

k  ((I.EQ.l). AID. (J.EQ.  0))  .OR. 

k  ((I.EQ.l) .AID. (J.EQ.-l))  )  RETURI 

IIITALIZE  OUTPUTS  TO  ZERO 

DO  5  1=1,3 
DO  5  K=1 ,2 

5  SURFIIT(I,J,I,K)=(0. ,0.) 

START  THE  OUTER  IITEGRATIOI  LOOP 
DO  50  JJ=1 ,8 

V=(I8(JJ)+2.D0*J+1.D0) *DELY/2 . DO 

IIITIALIZE  THE  IIIER  LOOP  STORAGE  MATRIX  ELEMEITS  TO  ZERO 

DO  10  1=1,3 
DO  10  K=l,2 

10  TERM(I,K)=(0. ,0.) 

START  IIIER  IITEGRATIOI  LOOP 

DO  20  11=1,8 

DEFIIE  U  AID  HAIK_TERM 

U=(I8(II)+2.DO*I-1.DO)*DELX/2.DO 

HAIK_TERM=HAIKO(KO*DSQRT(U*U+V*V)) 

START  PHYSICAL  BASIS  FUICTIOI  LOOP  AID  DEFIIE  EXP.TERM 


DO  20  1=1,3 

EXP_TERM=EXP ( -C J*G ( I) *U) 

C 

C  START  K  LOOP  AID  DEFIIE  TRIG.TERM 
C 

DO  20  K=1 ,2 

IF  (I.EQ.l)  THE! 

TRIG_TERM=DSII(F(I)*V) 
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ELSE 

TRIG_TERM=DCOS(F(I)*V) 

EIDIF 

C 

C  ADD  UP  THE  WEIGHTED  FUICTIOI  EVALUATION 
C 

TERM (I ,K)=TERM(I,K)+W8(II)* 
ft  TRIG_TERM*EXP_TERM*HAIK_TERM 

C 

C  EID  PHYSICAL  BASIS  FUICTIOI  AID  K  LOOPS 
C 

20  COITIIUE 

C 

C  ADD  WEIGHTED  IIRER  IITEGRATIOI  MATRIX  TO  OUTER  IITEGRATIOH  MATRIX 
C 

DO  40  1=1,3 
DO  40  K=1 ,2 

40  SURFIHT(I,J,I,K)  =  SURFIIT(I,J,I,K)  +  TERM(I ,K)*W8( JJ) 

C 

C  EID  OUTER  IITEGRATIOI  LOOP 
C 

50  COITIIUE 
C 

C  PERFORM  FIIAL  GAUSSIAI  QUADRATURE  MULTIPLICATIOI 
C 

DO  60  1=1,3 
DO  60  K=1 ,2 

60  SURFIIT ( I , J , I . K) =SURFIIT ( I , J , I , K) %DELX*DELY/4 . 

RETURI 

EID 
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2. 1-4-2  Subroutine  CORNERS 


C ********************************************************************** 


C  ** 
C  »»  SUBROUTIIE  COMERS  ««  ** 
C  ** 
C  THIS  SUBROUTIIE  CALCULATES  THE  ELEMEITS  OF  *♦ 
C  SURFIITd,  J,  I,  K)  THAT  DIRECTLY  TOUCH  THE  SIIGULAR  *♦ 
C  POUT.  THESE  ELEMEITS  ARE  GIVEI  BY  (I,J)=(0,0),  (0,-1),  ** 
C  (1,0),  AID  (1,-1).  FOR  EACH  CELL,  THE  IITEGRATIOI  OCCURS  ** 
C  OVER  THE  IORMAL  RECTAIGULAR  CELL  MIIUS  A  QUARTER-CIRCULAR  ** 
C  AREA  CEITERED  AT  THE  COMER  OF  THE  CELL  THAT  COIICIDES  ** 
C  WITH  THE  SIIGULAR  ORIGII.  ** 
C  ** 


C************************************ ************ ********************** 
C  »»  CALLED  BY:  MATRIX.FILL 

C  »»  SUBROUTIIES  CALLED:  IOIE 

C  »»  FUICTIOIS  CALLED:  HAIKO 

C  »»  COMMON  BLOCKS:  ALL,  GAUSSI AN.QUADRATURE ,  SURFACE.IITEGRALS 
C 

C  »»  IITERIAL  VARIABLES  «« 

C  R  RADIUS  OF  CIRCLE  AROUID  SIIGULAR  POUT 

C  U,  V  =  COORDIIATES  OF  POUT  AT  WHICH  FUICTIOI  IS  EVALUATED 

C  DURING  GAUSSIAN  QUADRATURE 

C  X  LIMIT  OF  IITEGRATIOI  DEFINED  BY  CIRCLE  AROUND  SINGULAR 

C  POINT 

C  SPAN  =  LENGTH  OF  IITEGRATIOI  INTERVAL  IN  U-DIMENSION 

C  TRIG.TERM  =  PART  OF  FUICTIOI  INTEGRATED 

C  HANK.TERM  =  PART  OF  FUICTIOI  INTEGRATED 

C  EXP.TERM  =  PART  OF  FUICTIOI  INTEGRATED 

C  TERM  =  TEMPORARY  STORAGE  MATRIX  OF  IIIER  INTEGRATIONS 

C  I,  J  =  DEFINE  COORDINATE  OF  LOWER-RIGHT  COMER  OF 

C  IITEGRATIOI  CELL 

C  I  PHYSICAL  BASIS  FUICTIOI  INDEX 

C  K  INTEGRATION  FUICTIOI  INDEX,  DEFINES  TRIG.TERM 

C  II  =  IIIER  GAUSSIAR  QUADRATURE  IIDEX 

C  JJ  =  OUTER  GAUSSIAI  QUADRATURE  IIDEX 

C 

C  »»  DATA  OUTPUT  TO  CALLING  ROUTIIE  «« 

C  SURFIITd, J.I.K)  WHERE  (I,J)  =  (0,0),  (0.-1),  (1.0),  OR  (1,-1) 

C 

C* ********************************************************************* 

SUBROUTIIE  COMERS 
C 

IMPLICIT  IOIE 

REAL*8  R,  U,  V,  X,  SPAI,  TRIG.TERM 

C0MPLEX*16  HAIK.TERM,  HAIKO,  EXP.TERM,  TERM(0: 1 ,-l : 0 , 1 : 3, 1 :2) 
INTEGER  I,  J.  I,  K,  II,  JJ 
C 

C  —  COMMON  BLOCKS - 

C 
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COMPLEX* 16  CJ 

REAL*8  PI,  EO,  HUO, ETA, KO, THETA, ER,H,V,F(3) ,G(3) , DELX, DELY 
ISTEGER  XIODES,  YVODES 

COMMOl  /  ALL  /  PI,  CJ,  EO,  HUO,  ETA,  KO,  THETA,  ER,  H,  W,  F,  G, 

*  XIODES,  YVODES,  DELX,  DELY 

C 

REAL*8  W8(l :8) ,  18(1:8),  W4(l:4),  14(1:4) 

COMMOl  /  GAUSSIAI.qUADRATURE  /  W8,  18,  v4,  n4 
C 

COMPLEX* 16  SURFIIT(-49:50,-5 :5, 1:3, 1:2),  SIIGIIT(4,3) 

COMMOl  /  SURFA CE_ IITEGRALS  /  SURFIIT,  SIIGIIT 
C 

C  DEFIIE  RADIUS  OF  CIRCLE  ABOUT  SIIGULAR  POUT.  DELX  >  DELY. 

C 

R=MII(DELY,  3./(DSQRT(ER)*K0)) 

IF  (R.GT.DELX)  THEN 

WRITE  (6,*)  'Error  in  CORKERS.  Slab  must  be  re-partitioned’ 
RETURN 
ENDIF 
C 

C  INITIALIZE  OUTPUT  TO  ZERO 
DO  10  1=0,1 
DO  10  J=-1,0 
DO  10  N=1 ,3 
DO  10  K=1 ,2 

10  SURFINT(I , J ,N,K)=(0. ,0. ) 

C 

C  OUTER  GAUSSIAN  QUADRATURE  INTEGRATION  LOOP 
DO  60  JJ=1 ,8 

V= (N8 ( J J ) + 1 ) *DELY/2 . DO 
X=O.DO 

IF  (V.LT.R)  X=DSQRT(R*R-V*V) 

SPAN=DELX-X 

C 

C  INITIALIZE  INNER  LOOP  STORAGE  MATRIX  TO  ZERO 
C 

DO  20  1=0,1 
DO  20  J=-1,0 
DO  20  N=1 ,3 
DO  20  K=1 ,2 

20  TERM(I, J,N,K)=(0. ,0.) 

C 

C  INVER  GAUSSIAN  QUADRATURE  INTEGRATION  LOOP 
C 

DO  40  11=1,8 

U=I8(II)*SPAN/2.  +  (DELX+X)/2. 
HANK_TERM=HAIKO(KO*DSQRT(U*U+V*V) ) 

C 

C  START  (I,J,N,K)  LOOP,  PERFORMING  FUNCTION  EVALUATION  FOR  EACH  CELL, 
C  PHYSICAL  BASIS  FUNCTION,  AND  TRIG.TERM. 

C 
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DO  30  1=0,1 
DO  30  J=-1,0 
DO  30  1=1,3 
DO  30  K=1 ,2 
C 

C  LOGIC  TO  DEFIHE  TRIG.TERM  AID  EXPJTERM 
C 

IF  (K.EQ.2)  THE! 

TRIG_TERM=DCOS(F(N)*V) 

ELSEIF  (J.EQ.-l)  THEM 

TRIG_TERM=DSIi(F(N)*(-V) ) 

ELSE 

TRIG_TERM=DSIN(F(N)*V) 

EIDIF 

C 

IF  (I.EQ.O)  THEN 

EXP_TERM=EXP (  CJ*G(N)*U) 

ELSE 

EXP_TERM=EXP(-CJ*G(N)*U) 

EHDIF 

C 

C  ADD  UP  THE  WEIGHTED  FUNCTION  EVALUATIONS 
C 

TERM (I , J,N,K)=TERM(I , J ,H,K)+W8(II)* 

*  TRIG_TERM*EXP_TERM*HANK_TERM 

C 

C  END  (I , J , I ,K)  LOOP 
C 

30  CONTINUE 

C 

C  END  INNER  GAUSSIAN  QUADRATURE  INTEGRATION  LOOP 
C 

40  CONTINUE 
C 

C  ADD  THE  WEIGHTED  INNER  INTEGRATION  TO  THE  OUTER  INTEGRATION  MATRIX 
C 

DO  SO  1=0,1 
DO  SO  J=-1,0 
DO  50  1=1 ,3 
DO  50  K=1 ,2 

SO  SURFIITd ,  J  ,N  ,K)  = 

t  SURFINT(I , J,N,K)+TERM(I , J,N,K)*W8(JJ)*SPAN/2. 

C 

C  END  OUTER  INTEGRATION  LOOP 
C 

60  CONTI  NUF. 

C 

C  PERFORM  THE  FINAL  GAUSSIAN  QUADRATURE  MULTIPLICATION 
C 

DO  70  1=0,1  . 

DO  70  J=-1,0 
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DO  70  1=1,3 
DO  70  K=1 ,2 

70  SURFIIT(I, J,I,K)=SURFIIT(I1 J 

RETURN 

END 


I,K)*DELY/2. 
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2.1.5  Subroutine  SINGULARJNTEGRAL 


C********************************************************************** 


c  ** 

C  »»  SUBROUTINE  SINGULAR.INTEGRAL  ««  ** 
C  ** 
C  THIS  SUBROUTINE  CALCULATES  THE  VALUE  OF  THE  SURFACE  ** 
C  AID  LIME  IITEGRALS  IV  THE  NEIGHBORHOOD  OF  THE  SIIGULAR  ** 
C  MATCH  POUT.  SIIGIIT  IS  ORGAIIZED  AS  SIIGIIT(K,I)  WHERE  ** 
C  ** 
C  K  =  1,  SURFACE  IITEGRAL  OVER  FULL  CIRCLE  *♦ 
C  2,  SURFACE  IITEGRAL  OVER  SEMICIRCLE,  EXPONENTIAL  ** 
C  COITAISS  +F(N)  ** 
C  3,  SURFACE  IITEGRAL  OVER  SEMICIRCLE,  EXPONENTIAL  ** 
C  COITAIIS  -F(I)  ** 
C  4,  LIME  IITEGRAL  ALONG  CONDUCTOR  ** 
C  N  =  PHYSICAL  BASIS  FUNCTION  INDEX  *♦ 
C  ** 
C  ALGORITHM  REPLACES  INTEGRANDS  BY  SMALL-ARGUMENT  POLYNOMIAL  ** 
C  EXPANSIONS  AND  INTEGRATES  TERM-BY-TERM.  ** 
C  ** 


C* **************************************************** ***************** 

C  CALLED  BY:  DEFINE 

C  SUBROUTINES  CALLED:  IONE 

C  FUNCTIONS  CALLED:  IITEGRAL_COSIIE,  FACT 

C  COMMON  BLOCKS:  ALL,  HANKEL,  SURFACE.INTEGRALS 

C 

C  »»  INTERNAL  VARIABLES  «« 

C  I  INDEX 

C  I  =  PHYSICAL  BASIS  FUNCTION  INDEX 

C  K,  KK  =  INDICES  OF  OUTER  AND  INNER  DO-LOOPS 
C  R  RADIUS  OF  IITEGRATIOI  REGION 

C  KOR  =  KO*R 

C  K0R32  =  K02  DIVIDED  BY  3,  QUANTITY  SQUARED 

C  LOGTERM=  NATURAL  LOGARITHM  OF  KOR 

C  ANGLE  =  VECTOR  CONTAINING  ANGLES  FOR  EACH  PBF 

C  GR2  =  VECTOR  CONTAINING  G(N)*R,  QUANTITY  SQUARED 

C  RATIO  =  CONVERGENCE  TEST  FOR  CURRENT  ITERATION 

C  OLD.RATIO  =  CONVERGENCE  TEST  FOR  PREVIDU  ITERATION 
C  TOL  =  CONVERGENCE  THRESHOLD 

C  JERKOR  =  CJ*SQRT(ER)*KOR 

C  INNER_SUM  =  RUNNING  SUM  FOR  INNER  DO-LOOP 

C  TERM  =  RUNNING  SUM  OF  OUTER  LOOP 

C 

C  »»  DATA  OUTPUT  TO  CALLING  ROUTINE  «« 

C  SINGINT  =  VALUE  OF  SURFACE  AND  LINE  INTEGRALS,  PASSED  VIA 
C  COMMON  BLOCK  "SURFACE. INTEGRALS" 

C 

C« ********************************************* ************************ 

SUBROUTIIE  SINGULAR.IITEGRAL 
C 
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IMPLICIT  IOIE 
IRTEGER  I,  I,  K,  KK 

REAL*8  R,  KOR,  K0R32,  LOGTERM,  AIGLE(3) ,  GR2(3), 

It  RATIO,  OLDJIATIO,  TOL,  IRTEGRAL.COSIRE,  FACT 
COMPLEX* 16  TERM,  USER. SUM,  JERXOR 
C 

C  —  COMMOR  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL *8  PI,  EO,  MUO, ETA, KO, THETA, ER,H,W,F(3) ,G(3) ,DELX,DELY 
IRTEGER  XRODES,  YRODES 

COMMOR  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  V,  F,  G, 
ft  XRODES,  YRODES,  DELI,  DELY 

C 

REAL*8  J0(0:6) ,  Y0(0:6),  BETA(0:6),  MAG0(0:6),  PHASE0(0:6) 
COMPLEX* 16  ALPHA(0:6) 

REAL*8  J1 CO : 6)  ,  Yl(0:6),  MAG1(0:6),  PHASE1(0:6) 

COMMOR  /  HAHKEL  /  JO,  YO,  ALPHA,  BETA,  MAGO,  PHASEO, 

A  Jl,  Yl,  MAGI,  PHASE1 

C 

COMPLEX* 16  SURFIHT(-49 :50,-6:5 , 1 : 3 , 1 : 2) ,  SIHGIHT(4,3) 

COMMOR  /  SURF ACE_ IRTEGRALS  /  SURFIRT,  SIRGIRT 
C 
C 
C 

C  IRITIALIZE  SIRGIRT  TO  ZERO  ARD  DEFIRE  SOME  IRTERRAL  VARIABLES. 

C 

TOL=l . 0D-8 

R=MIR (DELY ,3 . / (KO*DSQRT(ER) ) ) 

DO  10  H=1 ,3 

ARGLE(R)=DATAR2(F(R) ,G(R)) 

DO  10  K=1 ,4 

10  SIRGIRT(K,R)=(0. ,0.) 

K0R=K0*R 

LOGTERM=DLOG (KOR) 

JERKOR=CJ*DSQRT ( ER) *KOR 
K0R32=K0R*K0R/9 . 

C 

C  CALCULATE  SIRGIRT (1,*)  AS  THE  DOUBLE  SUMMATIOR  GERERATED  BY 
C  IRTEGRATIRG  THE  PRODUCT  OF  THE  SMALL-ARGUMEHT  POLYROMIAL 
C  APPROXIMATIORS  FOR  BESSEL  ARD  HARKEL  FURCTIORS  OF  ORDER  ZERO. 

C  LUCKILY,  SIRGIRT(1 , 1)=SIRGIRT( 1 , 2)=SIRGIRT(1 ,3)  . 

C 

DO  30  K=0,6 
TERM=(0. ,0.) 

DO  20  KK=0 ,6 

IRRER_SUM=IRRER_SUM+(K0R32**KK/ (K+KK+1))* 
ft  (ALPHA(KK)-CJ*BETA(KK)*(L0GTERM-1 . / (2* (K+KK+1)) )) 

20  CORTIRUE 

SIRGIRT(1 , 1)=SIRGIRT(1 , 1)  +  IRRER_SUM*JO(K)*(KOR32*ER)**K 
30  CORTIRUE 
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SIIGIIT( 1 , 1)=SIIGI1T(1 , 1)*PI*R*R 
SI1GI1T(1 ,2)=SI*GI*T(1 , 1) 

SIIGIITd  ,3)=SI*GI*Td .  D 
C 

C  SIIGIIT(2:3,*)  ARE  SURFACE  IITEGRALS  OVER  SEMICIRCLES,  AID  ARE 
C  CALCULATED  AS  DOUBLE  SUMMATIOIS.  WE  HAVE  TO  TEST  THE  OUTER 
C  SUMMATIOI  FOR  COIVERGEICE.  START  THE  OUTER  LOOP. 

C 

I=-l 

40  1=1+1 

C  THE  IHIER  LOOP  DOESM’T  DEPEID  01  I,  AID  ITS  RESULT  APPLIES 
C  TO  BOTH  SIIGIIT(2,*)  AID  SIIGI*T(3,*) 

C 

IIIER_SUM=(0. ,0.) 

DO  50  K=0,6 

IIIER_SUM=IIIER_SUM+ 
ft  (  K0R32++K  /  (I+2+K+2)  )  * 

ft  (  ALPHA (K) -CJ*BETA(K) * (LOGTERM  -  1 ./(I+2+K+2))  ) 

50  CONTIIUE 
C 

C  APPLY  THE  IIIER  LOOP  RESULT  TO  SIIGIIT(2:3,*) 

C 

DO  60  1=1,3 

TERM  =  IIIER.SUM  *  JERKOR++I  * 
ft  IITEGRAL_COSIIE(I , AIGLE (I) -PI , AIGLE (I) )/FACT(I ) 

SIIGIIT(2 ,1)  =  SIIGIIT(2,I)  +  TERM 
TERM  =  IIIER.SUM  *  JERKOR**!  * 
ft  IITEGRAL_COSIIE(I,-AIGLE(I)-PI,-AIGLE(I))/FACT(I) 

SIIGIIT(3,I)  =  SIIGIIT(3.I)  +  TERM 
60  COITIIUE 
C 

C  DUMP  OLD  RATIO  IITO  OLD.RATIO  AID  CALCULATE  IEW  RATIO;  GET  READY 
C  FOR  COIVERGAICE  TEST.  ASSUME  THESE  TERMS  COIVERGE  AT  THE  SAME 
C  RATE,  WHICH  IS  DOMIIATED  BY  THE  FACTORIAL. 

C 

OLD_RATIO=RATIO 

RATIO=ABS (TERM/SIIGIIT(2 , 1) ) 

C 

C  IF  COIVERGEICE  IOT  REACHED  II  OUTER  SERIES,  OR  LESS  THAI  FIVE 
C  TERMS  TAKER ,  GO  BACK  AID  FIID  AIOTHER  TERM. 

C 

IF  ((I.LE.5) .OR.(RATIO+OLD_RATIO.GT.TOL))  GOTO  40 
C 

C  DOIE  WITH  OUTER  SERIES;  MULTIPLY  BY  R**2 
C 

DO  70  *=1,3 

SIIGIIT(2,I)  =  SIIGIIT(2 ,1)  *  R  *  R 
SIIGIRT(3 ,1)  =  SI*GIIT(3 ,*)  *  R  ♦  R 
70  COITIIUE 
C 
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C  DOME  CALCULATIIG  SIIGIIT(2:3,*) .  01  TO  SIIGIIT(4,*) ,  THE  LIRE 

C  IITEGRAL.  R  CAI  IOW  BE  AS  LARGE  AS  3/KO,  IVFLUE1CIIG  SOME  OTHER 
C  IHTERVAL  VARIABLES. 

C 

R=MII ( 3 . /KO , DELX) 

KOR=KO*R 

LOGTERM=DLOG (KOR) 

K0R32=K0R*K0R/ 9 . 

DO  80  1=1,3 

80  GR2(I)  =  G(l)  *  G(l)  *  R  *  R 
IF  (GR2(l).EQ.O)  GR2(l)=l .D-10 
C 

C  START  THE  OUTER  LOOP;  AGAII  WE  MUST  TEST  FOR  COIVERGEICE. 

C 

I=-l 

90  1=1+1 
C 

C  START  THE  IIHER  LOOP 
C 

IIHER_SUM=(0. ,0.) 

DO  100  K=0,6 

IIIER_SUM=IIHER_SUM+ 

A  (  K0R32++K  /  (2+I+2+K+1)  )  * 

A  (  ALPHA (K)-CJ*BETA(K)*(LOGTERM  -  1 ./(2*I+2*K+1))  ) 

100  COHTIIUE 
C 

C  APPLY  THE  II1ER  LOOP  RESULT  TO  THE  CURREIT  OUTER  LOOP  TERM. 

C 

DO  110  1=1,3 

TERM=IIIER_SUM*(-GR2(I))+*I/FACT(2*I) 
SIIGIIT(4,I)=SI1GIIT(4,1)+TERM 
110  COHTIIUE 
C 

C  TEST  FOR  COIVERGEICE,  ASSUMIIG  ALL  SIIGIHT(4,*)  COIVERGE  AT 
C  THE  SAME  RATE. 

C 

RATIO=ABS ( TERM/SIIGIIT (4 , 1 ) ) 

IF  ((I.LE.B) .OR. (RATIO. GT.TOL))  GOTO  90 
C 

C  DOME  WITH  OUTER  SERIES;  MULTIPLY  BY  2R 
C 

DO  120  1=1,3 

SIIGIIT(4,I)=SIIGIIT(4,I)*2*R 
120  COHTIIUE 
C 

RETURI 

EID 
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2.2  Subroutine  IMPEDANCE-MATRIX 


C********************************************************************** 

c  ** 

C  »»  SUBROUTIIE  IMPED AICE_MATRIX  ««  ** 

C  ** 

C  THIS  SUBROUTIIE  GEIERATES  THE  IMPEDAICE  AID  VOLTAGE  ** 

C  MATRICES,  BOTH  UIIORMALIZED  AID  IORMALIZED.  UIIORMALIZED  ** 

C  MEAIS  IOM-SqUARE,  WHILE  IORMALIZATIOI  IIVOLVES  MATRIX  ** 

C  MULTIPLICATIOI  BY  THE  COIJUGATE  TRAISPOSE  OF  THE  ** 

C  UIIORMALIZED  IMPEDAICE  MATRIX.  ** 

C  ** 

C********************************************************************** 
C  CALLED  BY:  MAI I 

C  SUBROUTIHES  CALLED:  IORMALIZE,  IIVERT 

C  FUHCTIOHS  CALLED:  ZMI,  ARGD 

C  COMMOI  BLOCKS:  ALL,  IMPEDAICE,  MATCH.POIIT 

C 

C  »»  UTERI AL  VARIABLES  «« 

C  I  =  PHYSICAL  BASIS  FUHCTIOI  IIDEX 

C  K  =  IMPLIED  DO-LOOP  IIDEX 

C  POUT  =  VECTOR  COITAIIIIG  THE  MATCH  POUTS 

C  SHE  =  SIIE  OF  THETA,  THE  AIGLE  OF  IICIDEICE 

C  COSIIE=  COSIIE  OF  THETA 

C  EXP.TERM  =  PHASE  OF  IICIDEIT  WAVE  AT  MATCH  POIIT 

C  ZMIX  =  3X1  VECTOR  COITAIIIIG  UIIORMALIZED  IMPEDAICE  MATRIX 

C  ELEMEITS  FOR  THE  CURREIT  MATCH  POUT 

C 

C  »»  DATA  OBTAIIED  FROM  SUBROUTIIE  IORMALIZE  «« 

C  ZI  =  IORMALIZED,  3-by-3  IMPEDAICE  MATRIX 

C  VI  =  IORMALIZED.  3-by-l  VOLTAGE  MATRIX 

C 

C  »»  DATA  OBTAIIED  FROM  SUBROUTIIE  IIVERT  «« 

C  CUR=  3X1  CURREIT  VECTOR  (PHYSICAL  BASIS  FUICTIOI 

C  AMPLITUDES,  THE  FIIAL  FRUIT) 

C 

C* ********************************************************************* 

SUBROUTIIE  IMPEDAICE.MATRIX 
C 

IMPLICIT  IOIE 
IITEGER  I,  K 

REAL* 8  SHE,  COSIIE,  P0IIT(40) ,  ARGD 
COMPLEX* 16  EXP.TERM,  ZMIX (3) 

C 

C  —  COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MUO, ETA, KO, THETA, ER.H ,W,F(3) ,G(3) ,DELX ,DELY 
IITEGER  XIODES,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  W.  F,  G, 


27 


*  ZVOOES ,  YIODES,  DELI,  DELY 

C 

IITEGER  IUMMPS 
REAL* 8  Cl 

COMPLEX* 16  2(40,3),  V(40),  ZI(3,3),  VI(3),  CUR(3) 

COMMO*  /  IMPEDAICE  /  Z„  V  ,  ZI,  V*.  CUR,  Cl,  IUMMPS 
C 

IITEGER  I,  J 

REAL*8  X,  Y,  H02PY ,  H02PY2,  H02MY,  H02MY2, 
ft  V02PX,  W02PX2 ,  V02MX,  V02MX2 

COMMOI  /  MATCH.POIIT  /I,  J,  X,  Y,  H02PY,  H02PY2,  H02MY,  H02MY2, 
ft  V02PX ,  V02PX2,  W02MX,  W02MX2 

C 

C  IIITIALIZE  EVERYTHIIG  TO  ZERO 
C 

DATA  (V(II) ,11=1,40)  /  40* (0 . ,0 . )  / 

DATA  ((Z(II,JJ),II=1,40),JJ=1,3)  /  120*(0. ,0. )  / 

C 

C  DEFIIE  SHE  AID  COSIIE  OF  THETA. 

C 

SIIE=DSII (THETA) 

COSIIE=DCOS (THETA) 

C 

C  START  LOOP  TO  STEP  THROUGH  EACH  MATCH  POIMT 
C 

IUMMPS=0 

10  READ  (19,*,EID=30)  P0IIT(IUMMPS+1) 

C 

C  CALCULATE  MATCH  POIIT  X-  AID  Y-COORDIIATES  AID  IIDICES.  FOR  THIS 
C  IMPLEMEITATIOI,  J=0  BUT  WE  COULD  MODIFY  SO  THAT  MATCH  POIITS  WITHII 
C  THE  DIELECTRIC  COULD  BE  TAKEI . 

C 

I=IIT(POIIT (IUMMPS+1 ) * (XIODES+ 1 ) +0 . 5) 

IF  ( (I .LT. 1) . OR. (I .GT.XIODES) )  GOTO  10 
J=0 

IUMMPS=IUMMPS+ 1 
X=  -W/2.D0  +  I*DELX 
Y=  -H/2.D0  +  J*DELY 
W02PX=W/2.+X 
W02PX2=W02PX*W02PX 
W02MX=W/2.-X 
W02MX2=W02MX*W02MX 
H02PY=(H/2.+Y) 

H02PY2=H02PY*H02PY 

H02MY=(H/2.-Y) 

H02MY2=H02MY  *H02MY 
C 

C  CALCULATE  THE  IICIDEIT  FIELD  Ex  AT  THE  CURREIT  MATCH  POIIT. 

C 

EXP_TERM=EXP(-CJ*K0*(X*SIIE+H02MY*C0SIIE)) 

V(IUMMPS)=  COSIIE*EXP_TERM 
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c 

c 

c 

c 


c 

c 


c 

c 

c 

c 


c 


c 

c 

c 


CALL  ZMN  TO  GET  IMPEDANCE  MATRIX  ELEMENTS  FOR  CURRENT  MATCH  POINT. 
DUMP  RESULT  INTO  TEE  UNNORMALIZED  IMPEDANCE  MATRIX. 

CALL  ZMN  (ZMNX) 

DO  20  N=1 ,3 

Z(NUMMPS ,N)=ZMNX(N) 

20  CONTINUE 

GOTO  10 

30  WRITE  (20,40)  NUMMPS 

40  FORMAT  (IX,  12,  *  Match  Points:') 

WRITE  (20, ' (15(F4.2, IX)) ’ )  (POINT(K),  K=l, NUMMPS) 

NORMALIZE  IMPEDANCE  AND  VOLTAGE  MATRICES,  AND  SOLVE  THE  RESULTANT 
3X3  SYSTEM  OF  EQUATIONS. 

CALL  NORMALIZE 
CALL  INVERT  (3) 

WRITE  (20,50)  CN,  (ABS(CUR(N) ) ,  ARGD(CUR(N)) ,  N=l,3) 

50  FORMAT  (IX,  *  CN=\  F7.1,  2x,  3(»(».  F8.6,  F6.1,  ’)’)) 

DONE,  RETURN  TO  CALLING  ROUTINE 

RETURN 

END 
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S.S.l  Subroutine  ZMN 


O* ************************************************************ ******** 


c  ** 

C  »»  SUBROUTIIE  ZMl(ZMIX)  ««  ** 
C  ** 
C  THIS  FUNCTION  RETURHS  THE  IMPEDANCE  MATRIX  ELEMENTS  ** 
C  CORRESPOIDIIG  TO  THE  CURREIT  MATCH  POIVT.  ZMIX  CORRESPOIDS  ** 
C  TO  THE  Ex  IITEGRAL  EQUATIOI  ** 
C  ** 


C* ********************************************************************* 

C  CALLED  BY:  IMPED AICE.MATRIX 

C  SUBROUTIIES  CALLED:  SURF. I ITEGRALS ,  LIIE.IITEGRALS 

C  FUNCTIONS  CALLED:  COH JUGATE,  HAIKO,  HAHK1 
C  COMMO*  BLOCKS:  ALL,  FIELD.AMPLITUDES ,  MATCH.POIMT, 

C  SURFACE.IITEGRALS 

C 

C  »»  IITERIAL  VARIABLES  «« 

C  I  =  PHYSICAL  BASIS  FUICTIOI  IIDEX 
C  FOURJ  =  RECIPROCAL  OF  4*CJ 
C  EXP.TERM  =  COMPLEX  EXPOIEITIAL  IIVOLVIIG  X 

C  SHE  =  SIRE  IIVOLVIIG  Y 

C  COSINE  =  COSIIE  IIVOLVIIG  Y 

C  SIIE.H  =  SIRE  IIVOLVIIG  H 

C  COSIIE.H  =  COSIIE  IIVOLVIIG  H 

C  EXP1  *  COMPLEX  EXPOIEITIAL  IIVOLVIIG  V 

C  EXP2  =  COMPLEX  COIJUGATE  OF  EXPl 

C  PP  =  SqUARE  ROOT  OF  W02PX2  PLUS  H02PY2 

C  PM  =  SqUARE  ROOT  OF  V02PX2  PLUS  H02MY2 

C  MP  =  SqUARE  ROOT  OF  V02MX2  PLUS  H02PY2 

C  MM  =  SqUARE  ROOT  OF  W02MX2  PLUS  H02MY2 

C  HAIKO.PP  =  HARK EL  FUICTIOI  (ORDER  0)  IIVOLVIIG  PP 

C  HAIKO.PM  =  HAIKEL  FUICTIOI  (ORDER  0)  IIVOLVIIG  PM 

C  HAIKO.MP  =  HAIKEL  FUICTIOI  (ORDER  0)  IIVOLVIIG  MP 

C  HAIKO.MM  =  HAIKEL  FUICTIOI  (ORDER  0)  IIVOLVIIG  MM 

C  HAIKO.PP  =  HAIKEL  FUICTIOI  (ORDER  1)  IIVOLVIIG  PP 

C  HAIKO.MP  =  HAIKEL  FUICTIOI  (ORDER  1)  IIVOLVIIG  MP 

C 

C  »»  DATA  OBTAIIED  FROM  SUBROUTIIE  SURF.IITEGRALS  «« 

C  SURFACE  INTEGRALS  IMIA  AID  IMIB 

C 

C  »»  DATA  OBTAIIED  FROM  SUBROUTIIE  LIIE.IITEGRALS  «« 

C  LIIE  INTEGRALS  IMIG  THROUGH  IMIL 
C 

C********************************************************************** 
SUBROUTIIE  ZMI  (ZMIX) 

C 

IMPLICIT  IOIE 
IITEGER  I 

REAL*8  SIRE,  COSIIE,  SIIE.H,  COSIIE.H,  PP,  MP,  PM,  MM,  ER1 
COMPLEX* 16  FOURJ,  EXP.TERM,  EXPl,  EXP2,  HAIKO.PP,  HAIKO.PM, 
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ft  HAIKO.MP,  HAIKO.MM,  HAIK1.PP,  HAIK1.MP,  HAIKO,  HAIK1, 
ft  COI JUGATE,  SIIG.TERM 

COMPLEX* 16  ZMHX(3),  ZMIY(3),  ZMIZ(3),  IMHA(3),  IMNB(3), 
ft  IMIG(3) ,  IM1H(3) ,  IM1I(3),  IMIJ(3),  IMIX(3),  IMIL(3) , 

C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL* 8  PI,  EO,  MU0,ETA,K0,THETAtER,H,W,F(3) ,G(3) ,DELX,DELY 
I1TEGER  XIODES ,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  HUO,  ETA,  KO,  THETA,  ER,  H,  W,  F,  G, 
ft  XIODES,  YIODES,  DELX,  DELY 

C 

COMPLEX* 16  CX(3),  CY(3) ,  CZ(3) 

COMMOI  /  FIELD.AMPLITUDES  /  CX,  CY,  CZ 
C 

IITEGER  I,  J 

REAL*8  X,  Y,  H02PY,  H02PY2,  H02MY ,  H02MY2, 
ft  V02PX ,  V02PX2 ,  W02MX,  U02MX2 

COMMOI  /  MATCH.POIIT  /  I,  J,  X,  Y,  H02PY,  H02PY2,  H02MY,  H02MY2, 
ft  V02PX ,  V02PX2 ,  U02MX ,  V02MX2 

C 

COMPLEX* 16  SURFIIT(-49:50,-5:6, 1 :3, 1:2),  SIIGIIT(4,3) 

COMMOI  /  SURFACE_IITEGRALS  /  SURFIIT,  SIIGIIT 
C 
C 

C  DEFIIE  UTERI AL  VARIABLES  THAT  DO  IOT  DEPEID  01  I 
C 

F0URJ=1./(4.*CJ) 

ER1=ER-1. 

PP=DSQRT(W02PX2+H02PY2) 

PM=DSQRT(W02PX2+H02MY2) 

MP=DSQRT (W02MX2+H02PY2 ) 

MM=DSQRT(W02MX2+H02MY2) 

HAIKO_PP=HAIKO(KO*PP) 

HAIKO_PM=HAIKO (KO*PM ) 

HAIKO_MP=HAIKO(KO*MP) 

HAIKO_MM=HAIKO(KO*MM) 

HAIK1_PP=HAIK1 (KO*PP) 

HAIK1_MP=HAIK 1 (KO*MP) 

C 

C  CALL  SURF.IITEGRALS  AID  LIIE.IITEGRALS  TO  GET  TERMS  IMIA  THROUGH 
C  IMHI.  SURF_ IITEGRALS  GIVES  US  THE  SURFACE  IITEGRALS  OVER  THE  CROSS- 
C  SECTIOIAL  SURFACE  OF  THE  SLAB.  LIIE.IITEGRALS  GIVES  US  LIIE 
C  IITEGRALS  ACROSS  THE  FOUR  EDGES  OF  THE  SLAB. 

C 

CALL  SURF.IITEGRALS  (IMIA,  IMIB) 

CALL  LIIE.IITEGRALS  (IMIG,  IMIH,  IMII,  IMIJ,  IMIK,  IMIL) 

C 

C  START  LOOP  II  I  AID  DEFIIE  SOME  IITERJAL  VARIABLES  THAT  DEPEID  01  I 
C 
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DO  10  1=1,3 

EXP_TERH=EXP(-CJ*G(1I)*X) 

SIIE=DSII(F(I)* (Y+H/2 . ) ) 

C0SIIE=DC0S(F(I)*(Y+H/2. )) 

SIIE_H=DSII(F(I)*H) 

COSIIE_H=DCOS(F(I)*H) 

EXP1=EXP(  CJ*G(I)*W/2. ) 

EXP2=C0I JUGATE (EXP1 ) 

C 

C  CALCULATE  ZMIX(I),  USIIG  THE  CORRECT  SIIGULAR  TERN 
C 

SIIG_TERH=SIIE*SIIGIIT (1,1) 

IF  (J.EQ.O)  SIIG_TERH=(SIIGIIT(2,*)-SIIGI*T(3,I))/(2*CJ) 
ZMIX(I)=CX(I)*SI*E*EXP_TERM 

ZMIX(I)=ZMIX(I)-FOURJ*ER1*EXP_TERM*CX(I)*KO*KO* 

*  (COSIIE*IMHA(H)+SIIE*IMHB(H)+SIHG_TERM) 

ZHIX(I)=ZMIX(I)+0.25*ER1*EXP_TERM*CY(I)*G(I)* 
k  ( IMHK (I) -COSIIE_H* IMHL ( I ) ) 

ZMHX(I)=ZMIX(I)-F0URJ*ER1*CY(H)*((EXP1*HAIK0_PP 
k  -EXP2*HAHKO_MP)-COSIME_H*(EXPi*HAIKO_PM-EXP2*HAHKO_HM)) 

ZHIX ( I ) =ZMIX (H ) +FOUR J *ER1 *CX ( ■ ) *K0  * 
k  (  SIKE*(EXP1*W02PX*IMIG(I)+EXP2*W02HX*IMHH(M)) 

k  +C0SIIE*(EXP1*W02PX*IHII(I)+EXP2*V02MX*IM!IJ(I))  ) 

ZMK X (H) =ZHIX (I) +0 . 25*CZ(I) * (ETA/KO) *EXP_TERM* 
k  (K0*K0-G(I)*G(I))*IMIK(I) 

ZMHX(H)=ZMHX(I)+rGJRJ*CZ(I)*G(I)*(ETA/KO)* 
k  (EXPl*:.../1  j_PP-EXP2*HAIK0_KP) 

ZMHX(H)=',.’iPv.  A)-0.25*CZ(U)*ETA*(  EXP1*V02PX*HANK1  JPP/PP 

k  +  EXP2*W02MX *HAIK1_MP/MP  ) 

10  COITIIUL 
C 

RETURI 

PHD 
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2. 2. 1.1  Subroutine  SURF-INTEGRALS 


O ********************************************************************* 


c  ** 

C  »»  SUBROUTIIE  SURF.IITEGRALS  (IMIA,  IMIB)  ««  ** 
C  ** 
C  THIS  SUBROUTIIE  RETURIS  THE  THE  SURFACE  IITEGRALS  OVER  THE  ** 
C  DIELECTRIC  CROSS-SECTIOIAL  SURFACE  LESS  THE  CIRCULAR  ** 
C  REGIOI  ABOUT  THE  SIGULAR  MATCH  POUT.  IM1A  AID  IMIB  ** 
C  ARE  SIMPLY  THE  SUMS  OF  ELEMEITS  OF  MATRIX  ** 
C  SURFIIT,  WHICH  HAS  BEEI  PREVIOUSLY  CALCULATED.  ** 
C  ** 


C* ************************************************ ********************* 

C  CALLED  BY:  ZMH 

C  SUBROUTIVES  CALLED:  I0IE 

C  FUICTIOIS  CALLED:  IOHE 

C  COMMOI  BLOCKS:  ALL,  SURFACE.IITEGRALS,  MATCH.POIIT 
C 

C  »»  IITERIAL  VARIABLES  «« 

C  MM  =  OUTER  DO-LOOP  IIDEX 

C  II  =  IIHER  DO-LOOP  IIDEX 

C  I  =  PHYSICAL  BASIS  FUICTIOI  IIDEX 

C 

C  »»  DATA  OUTPUT  TO  CALLIVG  ROUTIIE  «« 

C  IMIA,  IMIB  =  3X1  VECTORS  COITAIIIMG  VALUES  OF  SURFACE  IITEGRALS 

C 

C*** *********************************************************** ******** 
SUBROUTIIE  SURF.IITEGRALS  (IMIA,  IMIB) 

C 

IMPLICIT  IOHE 

IITEGER  MM,  II,  I 

COMPLEX* 16  IMIA(3),  IMIB(3) 

C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MUO, ETA, KO, THETA, ER,H,W,F(3) ,G(3) .DELX.DELY 
IITEGER  XIODES,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  V,  F,  G, 
ft  XIODES,  YIODES,  DELX,  DELY 

C 

COMPLEX* 16  SURFIIT(-49:50,-5 :6, 1:3, 1:2),  SIIGIIT(4,3) 

COMMOI  /  SURFACE. IITEGRALS  /  SURFIIT,  SIIGIIT 
C 

IITEGER  I,  J 

REAL* 8  X,  Y,  H02PY,  H02PY2,  H02MY,  H02MY2, 
ft  V02PX,  V02PX2 ,  V02MX ,  V02MX2 

COMMOI  /  MATCH.POIIT  /  I,  J,  X,  Y,  H02PY,  H02PY2,  H02MY ,  H02MY2, 
ft  V02PX ,  V02PX2 ,  V02MX ,  V02MX2 

C 
C 
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C  INITIALIZE  OUTPUT  AID  RUNNING  SUMS  TO  ZERO 
C 

DO  10  1=1,3 

IMNA(N)=(0. ,0.) 

10  INNB(M)=(0. ,0.) 

C 

C  ADD  UP  SURFINT  ELEMENTS  TO  GENERATE  IHNA  AND  IMHB. 

C  MIDDLE  LOOP  STEPS  ACROSS  SLAB  HORIZONTALLY,  WHILE  INNER  LOOP  STEPS 
C  UP  SLAB  VERTICALLY. 

C 

DO  20  N=1 ,3 

DO  20  MM=1-I , XNODES+l-I 
DO  20  NN=-J , YNODES-J-1 

IMNA(N)=IMNA(N)+SURFINT(MM,NN ,N , 1) 
IMNB(N)=IMNB(N)+SURFINT(MM,NN,N,2) 

20  CONTINUE 
C 

RETURN 

END 
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2. 2. 1.2  Subroutine  LINE-INTEGRALS 


C* ********************************************************************* 


C  ** 
C  »»  SUBROUTINE  LINE.INTEGRALS  (IMNG,  IMNH,  IMNI,  ««  ** 
C  »»  IMIJ ,  IMNK,  I  MIL)  ««  ** 
C  ** 
C  THIS  SUBROUTIHE  RETURHS  LIIE  IITEGRALS  ALOHG  THE  EDGES  ** 
C  OF  THE  SLAB.  THE  LIIE  IITEGRALS  ARE  FUICTIOIS  OF  THE  ** 
C  MATCH  POUT  AID  THE  PHYSICAL  BASIS  FUICTIOI.  THE  ** 
C  SUBROUTIHE  USES  8-POINT  GAUSSIAN  QUADRATURE  TO  NUMERICALLY  ** 
C  INTEGRATE  THE  INTEGRAND.  ** 
C  ** 


C********************************************************************** 

C  CALLED  BY:  ZMN 

C  SUBROUTINES  CALLED:  NONE 

C  FUNCTIONS  CALLED:  HANKO 

C  COMMON  BLOCKS:  ALL,  GAUSSIAN.qUADRATURE,  MATCH.PDINT 

C 

C  »»  INTERNAL  VARIABLES  «« 

C  K  =  OUTER  DO  LOOP  INDEX 

C  II  =  INNER  DO  LOOP  INDEX 

C  H  =  PHYSICAL  BASIS  FUNCTION  INDEX 

C  U,  V  =  GAUSSIAN  QUADRATURE  NODE 

C  SINE  =  SHE  TERM  IN  INTEGRANDS 

C  COSINE  =  COSINE  TERM  IN  INTEGRANDS 

C  POS.SQRT  =  SQUARE  ROOT  TERM  CONTAINING  W02PX2  OR  H02PY2 

C  NEG_SQRT  =  SQUARE  ROOT  TERM  CONTAINING  W02MX2  OR  H02MY2 

C  POS_HANKO=  HANKEL  FUNCTION  (ORDER  0)  CONTAINING  POS_SQRT 

C  NEG_HANKO=  HANKEL  FUNCTION  (ORDER  0)  CONTAINING  NEG_SQRT 

C  P0S_HANK1=  HANKEL  FUNCTION  (ORDER  1)  CONTAINING  POS_SQRT 

C  NEG_HANK1=  HANKEL  FUNCTION  (ORDER  1)  CONTAINING  NEG.SQRT 

C 

C  »»  DATA  FROM  CALLING  ROUTINE  «« 

C  N  =  PHYSICAL  BASIS  FUNCTION  INDEX 

C 

C  »»  DATA  OUTPUT  TO  CALLING  ROUTINE  «« 

C  IMNG  THRU  IMNL  =  3X1  VECTORS  CONTAINING  VALUES  OF  LINE  INTEGRALS 

C 

C* ********************************************************************* 

SUBROUTINE  LINE. INTEGRALS  (IMNG,  IMNH,  IMNI,  IMNJ,  IMNK ,  IMNL) 

C 

IMPLICIT  NONE 
INTEGER  N,  II,  K 

REAL*8  U,  V,  SINE,  COSINE,  POS.SQRT,  NEG.SQRT 
COMPLEX* 16  POS.HANKO,  NEG.HANKO,  P0S.HAHK1,  IEG.HANK1,  HANKO. 
t  HANK1,  EXP .TERM 

COMPLEX* 16  IMNG(3),  IMIH{3),  IMII(3),  IMNJ(3),  IMNK(3),  IMNL(3) 

C 

C  —  COMMON  BLOCKS - 

C 
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COMPLEX* 16  CJ 

REAL* 8  PI,  EO,  MUO, ETA, KO, THETA, ER,H,W,F(3) ,G(3) ,DELX,DELY 
IITEGER  XIODES,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  V,  F,  G, 

*  XIODES,  YIODES.  DELX,  DELY 

C 

REAL*8  18(1:8),  18(1:8),  W4(l:4),  14(1:4) 

COMMOI  /  GAUSSIAI.QUADRATURE  /  V8,  18,  *4,  n4 
C 

IITEGER  I,  J 

REAL*8  X,  Y,  H02PY,  H02PY2,  H02MY ,  H02MY2, 
ft  V02PX ,  V02PX2,  V02MX,  U02MX2 

COMMOI  /  MATCH_POIIT  /  I,  J,  X,  Y,  H02PY,  H02PY2,  H02MY,  H02MY2, 
ft  V02PX,  W02PX2,  U02MX ,  W02MX2 

C 

COMPLEX* 16  SURFIIT(-49 :50, -S : 5 , 1 :3 , 1 :2) ,  SIIGIIT(4,3) 

COMMOI  /  SURFACE.IITEGRALS  /  SURFIIT,  SIIGIIT 
C 
C 

C  INITIALIZE  OUTPUT  TO  ZERO  AID  DEFIIE  SOME  IITERIAL  VARIABLES 
C 

DO  5  1=1,3 

IMIG(I)  =  (0.  ,  0.) 

IMIH(I)=(0. ,  0.) 

IMII(I)=(0. ,  0.) 

IMIJ(I)=(0 . ,  0.) 

IMHK(I)=(0 . ,  0.) 

IMIL(I)=(0. ,  0.) 

5  COITIIUE 
C 

C  IITEGRATE  ALONG  VERTICAL  EDGES.  THE  OUTER  LOOP  STEPS  UP  THE  EDGES 
C  OF  THE  SLAB. 

C 

DO  10  K=-J , YNODES-J-1 
DO  10  11=1,8 
C 

C  DEFINE  GAUSSIAN  QUADRATURE  IODE,  AID  SOME  IITERIAL  VARIABLES. 

C 

V= (18 ( I I ) +2 . *K+ 1 ) *DELY/2 . DO 
POS_SQRT=DSQRT (W02PX2+V*V ) 

IEG_SQRT=DSQRT(W02MX2+V*V) 

POS_HAIKO=HAIKO(KO*POS_SQRT) 

■EG_HAIKO=HAIKO (KO*IEG_SQRT) 

P0S_HAIK1 =HAIK 1 (KO*POS_SQRT ) 

IEG_HAIK1=HAIK1(K0*IEG_SQRT) 

C 

C  ADD  UP  THE  WEIGHTED  INTEGRAND  EVALUATED  AT  THE  CURRENT  IODE 
C 

DO  10  1=1,3 

SIIE=DSII(F(I)*V) 

COSIIE=DCOS(F(I)*V) 
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IMG (I) =IMIG (I)  +V8 ( II )  *COSINE*POS_HANK  1  /  POS.SQRT 
IMNH(N)=IMNH(N)+V8(II)*C0SIIE+IEG_HAIK1  /  IEG_SQRT 
IMII(I)=IHII(I)+V8(II)*  SIIE*P0S_HAIK1  /  POS.SQRT 
IMNJ(N)=IMNJ(N)+W8(II)*  SIIE*IEG_HAIK1  /  IEG_SQRT 
10  COITIIUE 
C 

C  MULTIPLY  BY  D ELY/2 
C 

DO  20  *=1,3 

IHIG(I) =IM*G (*) *DELY/2 . DO 
IMIH(I)=IMIH(I)*DELY/2.D0 
IMII (I) =IM*I (I) *DELY/2 . DO 
IMI J (I) =IM* J (I) *DELY/2 . DO 
20  COITIIUE 
C 

C  IHTEGRATE  ALOIG  THE  HORIZOITAL  EDGES.  THE  OUTER  LOOP  STEPS  FROM 
C  LEFT  TO  RIGHT  OIE  CELL  AT  A  TIME. 

C 

DO  40  K=-I,  XIODES-I 
DO  40  11=1,8 
C 

C  DEFIIE  GAUSSIAN  QUADRATURE  IODE,  AID  SOME  IITERNAL  VARIABLES 
C 

U= ( *8 ( I I ) +2 . *K+ 1 ) +DELX/2 . DO 
IEG_HAIKO=HAIKO (KO*DSQRT (U+U+H02MY2 ) ) 
P0S_HAIK0=HAIK0(K0*DSQRT(U*U+HQ2PY2) ) 

C  ADD  UP  THE  WEIGHTED  FUICTIOI  EVALUATIOIS  FOR  EACH  I.  CHECK  TO  SEE 
C  IF  IMNK  IS  HEAR  THE  SIIGULAR  POUT. 

C 

DO  40  1=1,3 

EXP_TERM=EXP ( -C J*G (I) *U) 
IMIL(I)=IM*L(I)+W8(II)*EXP_TERM*IEG_HANK0 
IF  ((J.EQ.O) .AID. ((K.EQ.-l) .OR. (K.EQ.O)))  GOTO  40 
IMIK(*)=IMIK(*)+W8(II)*EXP_TERM*P0S_HAIK0 
40  COITIIUE 
C 

C  MULTIPLY  BY  DELX/2,  AID  RETURI  IF  SINGULARITY  IS  IOT  01  THE 
C  CONDUCTOR 
C 

DO  60  1=1,3 

IMIK(I)=IMIK(I)*DELX/2 . 

IMIL(I)=IM*L(I)*DELX/2 . 

50  COITIIUE 

IF  (J.IE.O)  RETURI 
C 

C  HANDLE  THE  IITEGRATIOIS  BEAR  THE  SIIGULARITY,  IF  ON  THE  COIDUCTOR. 
C  DO  AS  MUCH  IUMERICAL  IITEGRATIOI  OF  IMIK  AS  IECESSARY . 

C  IF  DELX  IS  SMALL  EIOUGH  (DELX  <=  3/KO)  THEN  IOIE  IEEDED. 

C  EXPLOIT  THE  SYMMETRY  OF  THE  IITEGRAID . 

C 


IF  (K0*DELX.GT.3)  THEM 
DO  60  11=1,8 

U=0 . 5* ( (DELX-3 . /KO) *18 ( II ) +DELX+3 . /KO) 

DO  60  1=1,3 

IHIK(I)=IHIK(I)+(DELX-3./K0)*V8(II)* 

*  DCOS(G(I)*U)*HAIKO(KO*U) 

60  COITIIUE 
EIDIF 
C 

C  ADD  II  THE  AIALYTIC  IITEGRATIOI  AROUID  THE  SIIGULARITY. 
C 

DO  70  1=1,3 

IKIK(I)=IMIK(I)+SIIGIIT(4,I) 

70  COITIIUE 
C 

RETURI 

END 
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2.2.2  Subroutine  NORMALIZE 


C********************************************************************** 


C  ** 
C  »»  SUBROUTISE  NORMALIZE  ««  ** 
C  ** 
C  THIS  SUBROUTINE  NORMALIZES  THE  UNNORMALIZED  IMPEDANCE  AND  ** 
C  VOLTAGE  MATRICES  BY  MULTIPLYING  BOTH  BY  THE  CONJUGATE  ** 
C  TRANSPOSE  OF  THE  IMPEDANCE  MATRIX.  WHILE  THE  UNNORMALIZED  ** 
C  IMPEDANCE  AND  VOLTAGE  MATRICES  ARE  M-by-3  AND  M-by-1,  ** 
C  RESPECTIVELY,  THE  NORMALIZED  IMPEDANCE  AND  VOLTAGE  MATRICES  ** 
C  ARE  3-by-3  AND  3-by-l,  RESPECTIVELY.  ** 
C  ** 


C** ************************************************** ****************** 

C  CALLED  BY:  IMPEDANCE.MATRIX 

C  SUBROUTINES  CALLED:  NONE 

C  FUNCTIONS  CALLED:  CONJUGATE 

C  COMMON  BLOCKS:  ALL,  IMPEDANCE,  MATCH.POINT 

C 

C  »»  INTERNAL  VARIABLES  «« 

C  II,  JJ,  M  =  DO-LOOP  INDICES 

C 

C******** ***************************************************** ********* 
SUBROUTINE  NORMALIZE 
C 

IMPLICIT  NONE 
INTEGER  II,  JJ,  M 
COMPLEX* 16  CONJUGATE 
C 

C - COMMON  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MUO, ETA. KO, THETA, ER,H,W,F(3) ,G(3) .DELX.DELY 
INTEGER  XNODES,  YNODES 

COMMON  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER.  H,  W,  F,  G, 

*  XNODES,  YNODES,  DELX,  DELY 

C 

INTEGER  NUMMPS 
REAL* 8  CN 

COMPLEX* 16  Z(40,3),  V(40) ,  ZN(3,3),  VN(3),  CUR(3) 

COMMON  /  IMPEDANCE  /  Z,  V  ,  ZN,  VN,  CUR,  CN,  NUMMPS 
C 

INTEGER  I,  J 

REAL*8  X,  Y,  H02PY,  H02PY2,  H02MY,  H02MY2, 
ft  W02PX,  W02PX2 ,  W02MX ,  V02MX2 

COMMON  /  MATCH.POINT  /  I,  J,  X,  Y,  H02PY,  H02PY2,  H02MY,  H02MY2, 
ft  W02PX ,  V02PX2 ,  W02MX ,  W02MX2 

C 

C  GENERATE  THE  NORMALIZED  VOLTAGE  MATRIX  BY  MULTIPLYING  THE 
C  UNNORMALIZED  VOLTAGE  MATRIX  BY  THE  CONJUGATE  TRANSPOSE  OF  THE 
C  UNNORMALIZED  IMPEDANCE  MATRIX 
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c 

DO  10  11=1,3 

VI(II)=(0.,  0.) 

DO  10  M=1,NUMMPS 

VN(II)=VN(II)+CONJUGATE(Z(M,II))*V(N) 

10  CONTINUE 
C 

C  GENERATE  THE  NORMALIZED  IMPEDANCE  MATRIX  BY  MULTIPLYING  THE 
C  UNNORMALIZED  IMPEDANCE  MATRIX  BY  ITS  CONJUGATE  TRANSPOSE 
C 

DO  20  11=1,3 
DO  20  JJ=1 ,3 

ZH(II, JJ)=(0. ,  0.) 

DO  20  M=1 .NUMMPS 

ZN(II , JJ)=ZN(II , JJ)+CONJUGATE(Z(M, II) )*Z(M, JJ) 

20  CONTINUE 
C 

RETURN 

END 
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2.2.S  Subroutine  INVERT 


C*********************************************************** 


C  »»  SUBROUT  I  IE  IIVERT  ««  * 
C  * 
C  THIS  SUBROUTIIE  USES  GAUSSIAI  ELIMIIATIOI  TO  CALCULATE  * 
C  THE  IIVERSE  OF  THE  SQUARE  MATRIX  ’BIIPUT'.  THE  SUB-  * 
C  ROUTIHE  FIRST  PUTS  THE  IIPUT  VECTOR  I*TO  A  SCRATCH  VEC-  * 
C  TOR  CALLED  IIPUT,  THEM  USES  ROW  OPERATIOIS  TO  CHAIGE  * 
C  THE  MATRIX  TO  UPPER  TRIAIGULAR  FORM  AID  STORES  THE  PER-  * 
C  MUTATIOH  FACTORS  II  THE  ZERO  ELEMEITS  OF  THE  UPPER  TRI-  * 
C  A1GULAR  MATRIX.  IEXT  THE  IIVERSE  IS  CALCULATED  BY  * 
C  SOLVIIG  THE  TRIAIGULARIZED  MATRIX  AGAIIST  EACH  PERMUTA-  * 
C  TED  COLUMI  OF  THE  IDEITITY  MATRIX  WITH  DIMEISIOI  * 
C  'IDIM' .  THE  IIVERSE  OF  THE  BIIPUT  MATRIX  IS  PLACED  II  * 
C  THE  MATRIX  IIVERSE.  * 
C  * 


C**** ******************************************* ************* 

C  »  CALLED  BY:  SUBROUTIIE  SOLVIVRS  OR  SOLVITER 
C 

C  »  SUBROUTIIES  CALLED:  ♦*  IOIE  ** 

C 

C  »»  IIPUT  VARIABLES  «« 

C  BIIPUT  =  SQUARE  MATRIX  TO  BE  IIVERTED 
C  IDIM  =  DIMEISIOI  OF  THE  IIPUT  MATRIX 
C 

C  »»  IITERIAL  CALCULATIOI  VARIABLES  «« 

C  CZERO  =  COMPLEX  VALUE,  0.0  +  jO.O 

C  DEL  =  COLUMI  MATRIX  OF  CHAIGE  II  IDF.ITITY  MATRIX  COLUMI 
C  ELEMEIT  VALUES  RESULTING  FROM  TRIAIGULARIZATIOI 

C  OF  THE  IIPUT  MATRIX 

C  IIPUT  =  MATRIX  USED  AS  SCRATCH  FOR  IIVERSIOI 

C  IPIVOT  =  COLUMI  VECTOR  OF  ROW  PIVOT  IIFORMATIOI.  THE 

C  VALUE  OF  EACH  ELEMEIT  REPRESEITS  THE  ROW  I UMBER 

C  THAT  WAS  PIVOTED,  AID  THE  ELEMEIT  IUMBER  REPRE- 

C  SEITS  THE  LOCATIOI  FOR  THE  PIVOTED  VALUES 

C  COLMAX  =  ABSOLUTE  VALUE  OF  THE  MAXIMUM  VALUE  II  A  COLUMI 
C  OF  THE  IIPUT  VECTOR 

C  q  =  VARIABLE  USED  DURIIG  BACKSOLVIIG  OPERATIOI,  HAS 

C  THE  VALUE  OF  SUM  OF  PRODUCTS  OF  ALL  KIOWI  X  VALUES 

C  TIMES  THEIR  CORRESPOIDIIG  RHS  ELEMEIT,  SPECIFICALLY 

C  X(J)  =  (RHS(J)  -  Q)/IIPUT(I,J) 

C  RHS  =  COLUMI  MATRIX  USED  TO  STORE  THE  ELEMEITS  OF  KIOWI 
C  VALUES  USED  TO  SOLVE  FOR  THE  IDEITITY  MATRIX  COLUMI 

C  TEMP  =  TEMPORARY  STORAGE  VARIABLE  FOR  MATRIX  ROW  PIVOTIIG 
C  OPERATIOIS  AID  TRIAIGULARIZATIOI  OPERATIOIS 

C 

C  »»  OUTPUT  VARIABLES  «« 

C  IIVERSE  =  IIVERSE  OF  IIPUT  MATRIX 
C 

C* ************************************* ********************** 
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c 

SUBROUTIIE  I  WERT  (IDIN) 

C 

IMPLICIT  ROHE 

IITEGER  I,  J.  K,  L,  IDIM,  IPIV0T(40) 

COMPLEX* 16  RS(40),  DEL(40) ,  Q,  CZERO,  TEMP, 

+  ZI1V(40,40) ,  RHS,  ALPHAS,  IIPUT(40,40) 

REAL *8  COLMAX,  SUM.ZI,  SUM.ZIIV,  IORM_ZI,  I0RM.ZIIV 
C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI.  EO,  MUO, ETA, KO, THETA, ER,H,W,F(3),G(3).DELX,DELY 
IITEGER  XRODES ,  YIODES 

COMMOI  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  V,  F,  G. 
A  XIODES,  YIODES,  DELX,  DELY 

C 

IITEGER  IUMMPS 
REAL*8  Cl 

COMPLEX* 16  Z(40,3),  V(40) .  ZI(3,3),  VI(3),  CUR(3) 

COMMOI  /  IMPED AICE  /  Z,  V  ,  ZI,  VI,  CUR,  Cl,  IUMMPS 
C 
C 

CZERO  =  DCMPLX(0. ,0.) 

C 

C  TRAISFER  THE  ORIGIIAL  MATRIX  IITO  THE  SCRATCH  MATRIX. 

C 

DO  90  I  =  1 , IDIM 
DO  90  J  =  1 .IDIM 
90  IIPUT(I.J)  =  ZI(I,J) 

C 

C  OUTERMOST  DO  LOOP  REPEATS  FOR  EACH  COLUMI  OF  THE  MATRIX 
C  AS  THE  MATRIX  IS  BEIIG  REDUCED  TO  LOWER  TRIAIGULAR  FORM. 

C 

DO  100,  J  =  1 , (IDIM-1) 

COLMAX  =  ABS(IIPUT( J , J)) 

IPIVOT(J)  =  J 
C 

C  SEARCH  FOR  PIVOT  ROW. 

C 

DO  110,  K  =  J+l,  IDIM 

IF(ABS(IIPUT(K , J) ) . GT . COLMAI ) THEI 
COLMAX  =  ABS(IIPUT(K,J)) 

IPIVOT(J)  =  K 
EID  IF 

110  COITIIUE 
C 

C  IF  HEEDED,  PIVOT  ROW  J  WITH  ROW  IPIVOT(J) . 

C 

IF(IPIVOT( J) . IE. J)THEI 
DO  120,  K  =  J.IDIM 
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TEMP  =  IIPUT(J,K) 

IIPUT(J.K)  *■  IIPUT ( IPIVOT ( J )  , K ) 

IIPUT(IPIVOT(J) ,K)  =  TEMP 
120  COVTIIUE 

EID  IF 
C 

C  WITH  THE  MIX  IBS  VALUE  OF  THE  COLUM*  01  TOP,  REPLACE 
C  THE  REMAIII1G  EL EMEUS  OF  THE  C0LUM1  WITH  THE  PERMUTATI01 
C  FACTOR  OF  THE  ROW. 

C 

DO  130,  K  =  (J+l) , IDIM 

TEMP  =  DCMPLX(-1. .0.)  ♦  (I*PUT(K, J)/I1PUT(J,J)) 
IIPUT CK,J)  =  TEMP 
DO  130,  L  =  (J+l), IDIM 

130  IIPUT(K.L)  =  (TEMP  *  IIPUT(J,L>)  +  IIPUT(K.L) 

C 

C  EID  DO  LOOP  THAT  PUTS  IIPUT  MATRIX  II  UPPER  TRIAIGULAR  FORM. 
C 

100  COITIIUE 
C 

C  IEXT  CALCULATE  IIVERSE  MATRIX  BY  SOLVIIG  THE  EQUATIOH 
C  Ax  =  b  FOR  b  =  EACH  COLUMI  OF  THE  IDEITITY  MATRIX.  AID 
C  x  =  EACH  COLUMI  OF  THE  IIVERSE  MATRIX. 

C 

DO  150  J  =  1 , IDIM 
C 

C  FORM  THE  J-th  COLUMI  OF  THE  IDEITITY  MATRIX. 

C 

DO  155  K=1,IDIM 

155  RS(K)  =  CZERO 
RS(J)  *  DCMPLX(1. ,0.) 

C 

C  PIVOT  ELEMEITS  OF  IDEITITY  MATRIX  COLUMI  II  SAME 
C  ORDER  AS  THE  IMPEDAICE  MATRIX  ROWS  WERE  PIVOTED 
C  DURIIG  REDUCTIOH  TO  TRIAIGULAR  FORM. 

C 

DO  156,  K=1 , (IDIM-1) 

IF(IPIVOT(K) .IE.K)THEI 
TEMP  =  RS(K) 

RS(K)  =  RS(IPIVOT(K)) 

RS(IPIVOT(K) )  =  TEMP 
EIDIF 

156  COITIIUE 
C 

C  PERMUTATE  THE  RHS  MATRIX  ACCORDIIG  TO  THE  MULTIPLICATIOI 
C  FACTORS  STORED  II  THE  ELEMEITS  OF  THE  IIPUT  MATRIX. 

C 

DO  160,  K  =  2,  IDIM 
DEL(K)  =  CZERO 
DO  165,  L  =  1,  (K-l) 

DEL(K)  =  DEL(K)  +  (RS(L)  *  IIPUT(K.L)) 
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165  COITIIUE 

RS(K)  =  DEL(K)  +  RS(K) 

160  COITIIUE 
C 

C  HOW  BACKSOLVE  TO  FIID  THE  ELEMEITS  OF  THE  IIVERSE  MATRIX 
C  J-th  COLUMH.  ALGORITHM  ADAPTED  FROM  FIGURE  2.1,  01  PAGE 
C  29  OF  THE  TEXT  "IUMERICAL  AIALYSIS"  BY  JOHISOI  AID  REISS. 
C 


+ 

175 

+ 

170 
C 

160  COITIIUE 
C 

C  FI HD  CURREHT  VECTOR  AS  PRODUCT  ZIIV*VI 
C 

DO  200  1=1 , IDIM 
CUR(I)=CZERO 
DO  200  J=1,IDIM 

200  CUR(I)=CUR(I)+ZIHV(I,J)*VI(J) 

C 

C  CALCULATE  COIDITIOI  HUMBER 
C 

DO  300  1=1, IDIM 
SUM_ZH  =0.0 
SUM.ZIIV  =0.0 
DO  290  J=1 ,IDIM 

SUM_ZI  =  SUM_ZI+ABS(ZI(I,J>) 
SUM.ZIHV  =  SUM_ZIHV+ABS(ZIIV(I , J)) 

290  COITIIUE 

IORM_ZI=MAX ( HORM_ZH ,  SUM.ZI ) 
■ORM_ZIIV=MAX(IORM_ZIIV,  SUM.ZIIV) 

300  COITIIUE 

CI=IORM_ZI*IORM_ZIIV 
WRITE  (6,*)  'COIDITIOI  IUMBER  IS  ’ ,CI 
C 

1000  RETURH 
EID 


ZIIV(IDIM.J)  =  RS(IDIM)/IIPUT(IDIM, IDIM) 

DO  170,  K  =  1 , (IDIM-1) 
q  =  CZERO 
DO  175,  L  =  l.K 

q  =  Q  +  IIPUT ( ( IDIM-K) , ( IDIM- (L- 1 ) ) )  * 
ZIIV((IDIM-(L-1)),J) 

COHTIIUE 

ZIIV ( IDIM-K, J)  =  (RS(IDIM-K)  -  q)  / 

IIPUT( IDIM-K, IDIM-K) 

COHTIIUE 
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2.3  Subroutine  FIND.CURRENTS 


C*************************************************************** ******* 


c  ** 

C  »»  SUBROUTIIE  FIND.CURRENTS  ««  *♦ 
C  ** 
C  THIS  SUBROUTIIE  CALCULATES  THE  EQUIVALEIT  CURREITS  ALOHG  ** 
C  THE  STRIP  AID  ALONG  THE  MIDDLE  OF  THE  SLAB.  ** 
C  ** 


C* ********************************************************************* 

C  CALLED  BY:  PBFSTRIP 

C  SUBROUTINES  CALLED:  IOIE 

C  FUICTIOIS  CALLED:  IOIE 

C  COMMOI  BLOCKS:  ALL,  FIELD. AMPLITUDES ,  IMPEDANCE 

C 

C  »»  INTERNAL  VARIABLES  «« 

C  M  =  DO-LOOP  INDEX 

C  I  =  PHYSICAL  BASIS  FUNCTION  INDEX 

C  X  =  X-COORDINATE  OF  CURRENT-SAMPLING  POINT 

C  Y  =  Y-COORDINATE  OF  CURRENT-SAMPLING  POINT 

C  JEqX  =  X-DIRECTED  POLARIZATION  CURRENT 

C  JEQY  =  Y-DIRECTED  POLARIZATION  CURRENT 

C  JCX  =  X-DIRECTED  CONDUCTION  CURRENT 

C 

C ********************************************************************** 
SUBROUTINE  FIND.CURRENTS 
C 

IMPLICIT  NONE 
REAL*8  X,  Y 
INTEGER  N,  M 

COMPLEX* 16  JEQX,  JEQY,  JCX 
C 

C - COMMON  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL* 8  PI,  EO,  MUO, ETA, KO, THETA, ER,H,V,F(3) ,G(3) .DELX.DELY 
INTEGER  XNODES,  YNODES 

COMMON  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  V,  F,  G, 
ft  XNODES,  YNODES,  DELX,  DELY 

C 

COMPLEX* 16  CX(3),  CY(3) ,  CZ(3) 

COMMON  /  FIELD. AMPLITUDES  /  CX,  CY,  CZ 
C 

INTEGER  NUMMPS 
REAL* 8  CN 

COMPLEX* 16  Z(40,3) ,  V(40),  ZN(3,3),  VN(3),  CUR(3) 

COMMON  /  IMPEDANCE  /  Z,  V  ,  ZN,  VN,  CUR,  CN,  NUMMPS 
C 

C  STEP  ACROSS  STRIP  AND  SLAB,  0.05  WAVELENGTH  AT  A  TIME 
C 
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DO  30  M=0,W*20 

X=-W/2.  +  0.05*M 
Y=0 

JEQX=(0. ,0.) 

JEqy=(o. ,0.) 

JCX=(0. ,0.) 

c 

C  ADD  UP  THE  CURREITS  FOR  EACH  PBF 
C 

DO  10  1=1,3 

JEQX= JEQX+CUR(I) *CX(I )*DSII(F(I) *H/2 . ) *EXP(-CJ*G (I) *X) 
JEQY=JEQY+CUR(I)*CY(I)*DC0S(F(I)*H/2 . )*EXP(-CJ*G(I)*X) 
JCX=JCX+CUR(I)*CZ(I)*EXP(-CJ*G(I)*X) 

10  COITIIUE 

JEQX=JEQX*CJ*(K0/ETA)*(ER-1. ) 

JEQY=JEQY*CJ*(K0/ETA)*(ER-1. ) 

C 

C  WRITE  RESULTS  TO  ASCII  AID  UIFORMATTED  FILES 
C 

WRITE  (31)  SIGL(H) ,  SIGL(ABS( JEQX)) , 
k  SBGL(ABS( JEQY) ) ,  SIGL(ABS( JCX) ) 

WRITE  (21,20)  X,  ABS(JEqX),  ABS(JEqY),  ABS(JCX) 

20  FORMAT  (IX,  F6.2,  3(31,  E20.10)) 

30  COITIIUE 
C 

RETURI 

EID 
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£.4  Subroutine  ERRORS 


C********************************************************************** 


C  ** 
C  »»  SUBROUTIIE  ERRORS  ««  *♦ 
C  ** 
C  THIS  SUBROUTIIE  CALCULATES  THE  RMS  AID  AVERAGE  RELATIVE  ** 
C  ERRORS  ASSOCIATED  WITH  THE  RESIDUAL  OF  THE  LEAST-SQUARES  ** 
C  IORMALIZATIOI.  ** 
C  ** 


C* ***************************************** **************************** 

C  CALLED  BY:  PBFSTRIP 

C  SUBROUTIIES  CALLED:  IOIE 

C  FUHCTIOHS  CALLED:  IOHE 

C  COMMOI  BLOCKS:  IMPED AICE 

C 

C  »»  IITERIAL  VARIABLES  «« 

C  RMS.ERROR  =  RMS  ERROR 

C  AVG_REL_ERROR  =  AVERAGE  RELATIVE  ERROR 

C  ERROR  =  IITERMEDIATE  TERM 

C  II  =  DO-LOOP  IIDEX 

C 

O********************************************************************* 

SUBROUTIIE  ERRORS 
C 

IMPLICIT  IOHE 

REAL*8  ERROR,  RMS_ERROR,  AVG_REL_ERROR 
IITEGER  II 
C 

C - COMMOI  BLOCKS - 

C 

IITEGER  IUMMPS 
REAL*8  Cl 

COMPLEX* 16  Z(40,3),  V(40) ,  ZI(3,3),  VI(3),  CUR(3) 

COMMOI  /  IMPEDAICE  /  Z,  V  ,  ZI,  VI,  CUR,  Cl,  IUMMPS 
C 

C  CALCULATE  THE  ROOT-MEAI-SQUARE  (RMS)  ERROR 
C 

RMS_ERR0R=0 . 0 
DO  10  11=1, IUMMPS 

ERR0R=ABS(Z(II,1)*CUR(1)+Z(II,2)*CUR(2)+ 

A  Z(II,3)*CUR(3)-V(II))**2 

RMS_ERROR=RMS_ERROR+ERROR 
10  COITIIUE 

RMS_ERROR=DSQRT (RMS.ERROR/IUMMPS ) 

C 

C  CALCULATE  THE  AVERAGE  RELATIVE  ERROR 
C 

AVG_REL_ERROR=0 . 0 
DO  20  11=1, IUMMPS 
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ERROR=ABS( (Z(II , 1)*CUR(1)+Z(II ,  2)*CUR(2)  + 

*  Z(II,3)*CUR(3)-V(II))/V(II)) 

AVG_REL_ERROR=AVG_REL_ERROR+ERROR 
20  COITIIUE 

AVG_REL_ERROR=AVG_REL_ERROR/*UMMPS 
WRITE  (20,30)  RMS.ERROR,  AVG_REL_ERROR 
30  FORMAT  (’  RMS  Error  =  F8.6,  »,  Avg  Rel  Error  =  F8. 

RETURI 

EID 
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2.5 


Subroutine  RADIATE 


C********************************************************************** 


C  ♦  * 
C  »»  SUBROUTIBE  RADIATE  «« 

C  ** 
C  THIS  SUBROUTIBE  CALCULATES  THE  RMS  ABD  AVERAGE  RELATIVE  ** 
C  ERRORS  ASSOCIATED  WITH  THE  RESIDUAL  OF  THE  LEAST-SQUARES  ** 
C  BORMALIZATIOB .  ** 
C  ** 


C*************** *************** ********* ******************** *********** 

C  CALLED  BY:  PBFSTRIP 

C  SUBROUTIBES  CALLED:  IOBE 

C  FUBCTIOBS  CALLED:  SIBC 

C  COMMOB  BLOCKS:  ALL,  FIELD.AMPLITUDES ,  IMPEDANCE 
C 

C  »»  INTERNAL  VARIABLES  «« 

C  B  =  PHYSICAL  BASIS  FUNCTION  INDEX 
C  ANGLE  =  BISTATIC  ANGLE  DO-LOOP  INDEX 

C  TERM1 ,  TERM2  =  INTERMEDIATE  TERMS 

C  TERM  =  3X1  VECTOR  CONTAINING  INTERMEDIATE  TERMS 

C  EXPTERM  =  INTERMEDIATE  TERM  INVOLVING  AN  EXPONENTIAL 

C  INTEGRAL  =  PROPORTIONAL  TO  FIELD  INTENSITY 

C  PHI  =  BISTATIC  ANGLE  IN  RADIANS 

C  SINE  *  SINE  OF  PHI 

C  COSINE  =  COSINE  OF  PHI 

C  SINEH  =  3X1  VECTOR  CONTAINING  INTERMEDIATE  TERMS 

C  COSINEH  =  3X1  VECTOR  CONTAINING  INTERMEDIATE  TERMS 

C  SIGMA  =  TOTAL  SCATTERING  WIDTH  IN  dB 

C  PBF.SIGMA  =  3X1  VECTOR  CONTAINING  SCATTERING  WIDTH 

C  CONTRIBUTED  BY  THE  INDIVIDUAL  PBFs 

C 

C* ********* ************************************************************ 

SUBROUTINE  RADIATE 
C 

IMPLICIT  NONE 
INTEGER  J.  N,  ANGLE 

REAL*8  PHI,  SINE,  COSINE,  SIHEH(3),  C0SINEH(3), 

*  SIGMA,  PBF_SIGMA(3) ,  SINC 

C0MPLEX*16  TERM1 ,  TERM2 ,  TERM ( 3 ) ,  EXPTERM,  INTEGRAL 
C 

C - COMMON  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MUO, ETA, KO, THETA, ER, H, W, F(3) ,G(3) .DELX.DELY 
INTEGER  XNODES,  YNODES 

COMMON  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  W,  F,  G, 
ft  XNODES,  YNODES,  DELX,  DELY 

C 

COMPLEX* 16  CX(3),  CY(3),  CZ(3) 
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COMMON  /  FIELD, AMPLITUDES  /  CX,  CY,  CZ 
C 

INTEGER  NUMMPS 
REAL*8  CN 

COMPLEX* 16  Z(40,3) ,  V(40),  ZN(3,3),  VN(3),  CUR(3) 

COMMON  /  IMPEDANCE  /  Z,  V  ,  ZN,  VN,  CUR,  CN,  NUMMPS 
C 

C  CALCULATE  SOME  INTERMEDIATE  TERMS 
C 

DO  10  N=1 ,3 

SINEH(N)=DSIN(F(N)*H) 

COSINEH(N)=DCOS(F(i)*H) 

10  CONTINUE 
C 

C  STEP  THROUGH  EACH  BISTATIC  ANGLE,  ONE  DEGREE  AT  A  TIME 
C 

DO  40  AHGLE=0,  180 
PHI=ANGLE*PI/ 180 . 

SINE=DSIN(PHI) 

COSI NE=DCOS ( PHI ) 

EXPTERM=EXP ( CJ*KO*H*SINE) 

INTEGRAL=(0 . ,0.) 

C 

C  CALCULATE  THE  CONTRIBUTION  FROM  EACH  PBF  AND  SUM  THEM 
C 

DO  20  N=1 ,3 

TERM1=EXPTERM*(CJ*K0*SINE*SINEH(N)-F(N)*C0SINEH(N))+F(N) 
TERM2=EXPTERM* (C J*KO*SINE*COSIHEH (N) +F(N ) *SINEH(N) ) 

*  -CJ*KO*SINE 

TERM(N)=CX(N)*SINE*TERM1  -  CY(N)*C0SINE*TERM2 
TERM(N)=CJ*(K0/ETA)*(ER-1 . )*TERM(N)/ 
ft  (F(N)*F(N)-KO*KO*SINE*SINE)+CZ(N)*SINE 

TERM(N)=TERM(N)*SINC((K0*C0SINE-G(N))*W/2.)*W 
TERM(N)=TERM(N)*CUR(N)*EXP(-CJ*K0*H*SINE/2. ) 
INTEGRAL=INTEGRAL+TERM(N) 

IF  (TERM(N) .EQ.O*CJ)  THEN 
PBF_SIGMA(N)=- 100.0 
ELSE 

PBF_SIGMA(N)=10*DL0G10( . 25*K0*ETA*ETA*ABS(TERM(N) )**2) 
ENDIF 
20  CONTINUE 

SIGMA= 10. *DLOG 10(0. 26*K0*ETA*ETA*ABS (INTEGRAL) **2) 

C 

C  WRITE  THE  RESULTS  TO  ASCII  AND  UNFORMATTED  FILES 
C 

WRITE  (32)  SNGL(ANGLE),  (SNGL(PBF_SIGMA(N)) ,  N=l,3), 
ft  SNGL(SIGMA) 

WRITE  (22,30)  ANGLE,  (PBF_SIGMA(N) ,  N=i,3).  SIGMA 
30  FORMAT  (IX,  14,  4(3X,  F10.4)) 

40  CONTINUE 
C 
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RETURI 

EID 
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III.  Auxiliary  Functions 

S.l  Function  HANKO 


C** ****************************************************** ************** 


C  ** 
C  »»  FUICTIOI  HAIKO(X)  ««  ** 
C  ** 
C  THIS  FUICTIOI  RETURNS  THE  HAIKEL  FUICTIOI  OF  THE  SECOID  ** 
C  KIID  OF  ORDER  0  OF  POSITIVE,  REAL  ARGUMEIT  X.  THE  ** 
C  ALGORITHM  IS  BASED  01  THE  SERIES  EXPAISIOIS  FOUID  ** 
C  II  AHS-55.  ** 
C  ** 


C* ********************************************* ************************ 
C  SUBROUTIIES  CALLED:  IOIE 

C  FUICTIOIS  CALLED:  IOIE 

C  COMMON  BLOCKS:  HANKEL 

C 

C  »»  IITERIAL  VARIABLES  «« 

C  XX  X*X/9  IF  X  <=  3,  OR  3/X  IF  X  >  3 

C  LOGTERM  =  IATURAL  LOAGRITHM  OF  X 

C  MAG  =  MAGIITUDE  OF  HAIK1  WHEN  X  >  3 
C  FHASE  =  PHASE  OF  HAIK1  WHEN  X  >  3 

C  I  DO-LOOP  IIDEX 

O******************* ************************************************** 

COMPLEX* 16  FUICTIOI  HAIKO(X) 

C 

IMPLICIT  IOIE 

REAL*8  X,  XX,  LOGTERM,  MAG,  PHASE 
IITEGER  I 
C 

C - COMMOI  BLOCKS - 

C 

REAL* 8  J0(0 : 6) ,  Y0(0:6),  BETA(0:6),  MAGO(0:6),  PHASE0(0:6) 
COMPLEX* 16  ALPHA(0:6) 

REAL* 8  Jl(0:6) ,  Yl(0:6),  MAG1(0:6),  PHASE1(0:6) 

COMMOI  /  HAIKEL  /  JO,  YO,  ALPHA,  BETA,  MAGO,  PHASEO, 

*  Jl,  Yl,  MAGI,  PHASE1 

C 
C 

C  X  MUST  BE  A  POSITIVE,  REAL  IUMBER.  RETURI  ZERO  IF  X=0. 

C 

HAIK0=(0.  ,  0.) 

IF  (X.LT.O.O)  X=-X 
IF  (X.Eq.O.O)  RETURI 
C 

C  IF  X  IS  LESS  THAI  OR  EqUAL  TO  3,  USE  THE  SMALL- ARGUMEIT  POLYIOMIAL 
C  EXPAISIOI  FOUID  II  AMS-55  EqUATIOIS  9.4.1  AID  9.4.2 
C 
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IF  (X.LE.3)  THEM 
XX=X*X/9.D0 
LOGTERM=DLOG(X) 

HANKO=ALPHA (6) - (0 . .  1 . ) ♦BETA (6 ) ♦LOGTERM 
DO  10  1=5, 0,-1 

HANKO=HANKO*XX  +  ALPHA(I)-(0 . ,  1 . )*BETA (I) ♦LOGTERM 
10  COKTIIUE 
C 

C  IF  X  IS  GREATER  THAI  3,  USE  THE  LARGE-ARGUMENT  ASYMPTOTIC 
C  EXPANSION  FOUID  II  AMS-55  EQUATIOI  9.4.3 
C 

ELSE 

XX=3.D0/X 

MAG=MAG0(6) 

PHASE=PHASE0(6) 

DO  30  1=5, 0,-1 

PHASE=PHASE*XX+PHASEO ( I ) 

MAG=M AG*XX+M AGO ( I ) 

30  COITIIUE 

HAHK0=EXP(-(0. ,  1 . )*(PHASE+X) )*MAG/DSQRT(X) 

EHDIF 

C 

RETURN 

END 
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S.2  Function  HANKl 


C********************************************************************** 


C  ** 
C  »»  FUICTIOI  HANKl (X)  ««  ** 
C  ** 
C  THIS  FUNCTION  RETURNS  THE  HAIKEL  FUICTIOI  OF  THE  SECOID  ** 
C  KIND  OF  ORDER  1  OF  POSITIVE,  REAL  ARGUMENT,  X.  THE  ** 
C  ALGORITHM  IS  BASED  01  THE  SERIES  EXPAISIOIS  FOUID  ** 
C  II  AMS-5S .  ** 
C  ** 


c* ************************************************************* ******** 

C  SUBROUTINES  CALLED:  IOIE 

C  FUICTIOIS  CALLED:  IOIE 

C  COMMON  BLOCKS:  ALL,  HAIKEL 

C 

C  »»  INTERNAL  VARIABLES  «« 

C  XX  =  X/3  IF  X  <=  3,  OR  3/X  IF  X  >  3 

C  REAL  =  REAL  PART  OF  HANKl  WHEN  X  <=  3 

C  IMAG  =  IMAGINARY  PART  OF  HANKl  VHEI  X  <=  3 

C  MAG  =  MAGNITUDE  OF  HANKl  WHEN  X  >  3 

C  PHASE  =  PHASE  OF  HANKl  VHEI  X  >  3 

C  I  DO-LOOP  INDEX 

C ********************************************************************** 
C0MPLEX*16  FUICTIOI  HANKl (X) 

C 

IMPLICIT  IOIE 

REAL*8  X,  XX.  REAL,  IMAG,  MAG,  PHASE 
INTEGER  I 
C 

C - COMMOI  BLOCKS - 

C 

COMPLEX* 16  CJ 

REAL*8  PI,  EO,  MU0,ETA,K0,THETA,ER,H,W,F(3) ,G(3) ,DELX,DELY 
IITEGER  XNODES,  YIODES 

COMMON  /  ALL  /  PI,  CJ,  EO,  MUO,  ETA,  KO,  THETA,  ER,  H,  V,  F,  G, 

*  XNODES,  YIODES,  DELX,  DELY 
C 

REAL*8  J0(0:6) ,  Y0(0:6),  BETA(0:6),  MAG0(0:6),  PHASE0(0:6) 
COMPLEX* 16  ALPHA (0:6) 

REAL*8  J 1 (0 : 6 ) ,  Yl(0:6),  MAG1(0:6),  PHASE1(0:6) 

COMMON  /  HAIKEL  /  JO,  YO,  ALPHA,  BETA,  MAGO,  PHASEO, 

*  Jl,  Yl,  MAGI,  PHASE1 
C 

C 

C  X  MUST  BE  A  POSITIVE,  REAL  NUMBER.  RETURN  ZERO  IF  X=0. 

C 

HANK1=(0 .  ,  0.) 

IF  (X.LT.O.O)  X=-X 
IF  (X.EQ.O.O)  RETURN 
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c 

C  IF  X  IS  LESS  THAI  3,  USE  THE  SMALL- ARGUMENT  SERIES  EXPANSION 
C  FOUND  IN  AMS-SS  EQUATIONS  9.4.4  AND  9.4.5 
C 

IF  (X.LE.3)  THEN 
XX=(X*X/9.D0) 

REAL=O.DO 

REAL=J1(6) 

IMAG=Y1(6) 

DO  10  1=5, 0,-1 

REAL=REAL*XX+J1(I) 

10  IMAG=IMAG*XX+Y 1(1) 

REAL=REAL*X 

IMAG= (2 . DO/PI) *DLOG (X/2 . DO)*REAL+IMAG/X 
HANK1=REAL-CJ*IMAG 
C 

C  IF  X  IS  GREATER  THAN  3,  USE  THE  LARGE- ARGUMENT  ASYMPTOTIC 
C  FORM  FOUND  IN  AMS-55  EQUATION  9.4.6 
C 

ELSE 

XX=3.D0/X 

MAG=MAG1(6) 

PHASE=PHASE1 (6) 

DO  30  1=5,0, -1 

PHASE=PHASE*XX+PHASE1 ( I ) 

30  MAG=MAG*XX+MAG1 (I) 

PHASE=PHASE+X 

HANK 1= (MAG/DSqRT(X) ) *EXP (-CJ*PHASE) 

ENDIF 

C 

RETURN 

END 
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5.5  Function  HANK2 


£*»***************************************•**************************** 


c  ** 

C  »»  FUICTIOI  HAVK2(X)  ««  ** 
C  M 
C  THIS  FUICTIOI  RETURIS  THE  HAVKEL  FUICTIOI  OF  THE  SECOID  ♦* 
C  KIID  OF  ORDER  TWO  OF  REAL  ARGUHEIT  X.  ** 
C  ** 


C********************************************************************** 

C  FUICTIOIS  CALLED:  HAIKO,  HAIK1 
C 

COMPLEX* 16  FUICTIOI  HAIK2(X) 

C 

REAL*8  X 

COMPLEX* 16  HAIKO,  HAIK1 
C 

C  X  MUST  BE  A  POSITIVE,  REAL  lUMBER.  RETURI  ZERO  IF  X=0. 

C 

HAIK2=(0. ,0.) 

IF  (X.LT.O.O)  X=-X 
IF  CX.Eq.O.O)  RETURI 
C 

C  USE  A  RECURREICE  RELATIOI  TO  GEIERATE  HAIK2 
C 

HAIK2=(2 . /X)*HAIK1 (X)-HAIKO(X) 

C 

RETURI 

EID 
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Function  BINOMIAL 


3.4 


C** ******************************************************************** 


c  ** 

C  »»  FUICTIOI  BIIOMIAL(I.J)  ««  ** 
C  ** 
C  THIS  FUICTIOI  RETURIS  THE  BIIOMIAL  COEFFICIEIT  OF  ** 
C  "I  CHOOSE  JM  WHERE  I  >=  J.  FORMULA  FROM  AMS-65.  ** 
C  ** 


C* ********************************************** *********************** 
REAL*8  FUICTIOI  BIIOMIAL(I , J) 

C 


IITEGER  I,  J 
C 

C  IF  I<]  THEM  RETURI  ZERO  AID  FLAG  THE  ERROR 
C 


IF  (I.LT.J)  THEM 

WRITE  (6.*)  'Illegal  arguments  in  FUICTIOI  BIIOMIAL’ 
BIIOMIAL=O.DO 
RETURI 
EIDIF 


C 

C  IF  J=0  OR  J=I,  RETURI  1 
C 

BII0MIAL=1 .DO 

IF  ((J.Eq.O).OR.(J.Eq.I))  RETURI 
C 

C  CALCULATE  BIIOMIAL,  USIIG  AMS-55'S  EqUATIOI  24. 1.1. C 
C 

DO  10  11=1, J 

10  BIIOMIAL=BIIOMIAL*DBLE (I-II+l ) /DBLE ( J-II+ 1 ) 

C 

RETURI 

EID 
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5.5  Function  FACT 


C********************************************************************** 


c  ** 

C  »»  FUNCTION  FACT(X)  ««  ** 

C  ** 

C  THIS  FUICTIOI  RETURNS  THE  FACTORIAL  OF  A  IOI-IEGATIVE  ** 

C  IITEGER,  X. 

C  ** 


C*** ************************************************ ******************* 

REAL*8  FUICTIOH  FACT(X) 

C 

IITEGER  X 
C 

C  FACT(0)=1 
C 

FACT=1 . 

IF  (X.EQ.O)  RETURN 
C 

C  CHECK  FOR  UPPER  LIMIT  OH  X 
C 

IF  (X.GT.30)  THEI 

WRITE  (6,*)  'Overfloa  in  SUBROUTINE  FACT' 

FACT=9 . 9999E+32 
RETURN 
ENDIF 
C 

C  C^ECK  FOR  ILLEGAL  NEGATIVE  X 
C 

IF  (X.LT.O)  THEN 

WRITE  (6,*)  'Negative  argument  in  SUBROUTINE  FACT' 

FACT=O.DO 

RETURN 

ENDIF 

C 

C  CALCULATE  FACT 
C 

DO  10  1=1, X 

10  F  ACT=FACT*DBLE ( I ) 

C 

RETURN 

END 
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3.6  Function  INTEGRAL-COSINE 


C********************************************************************** 


C  ** 
C  »»  FUICTIOI  IITEGRAL_COSIIE(I ,  A,  B)  ««  ** 
C  ** 
C  THIS  FUICTIOI  CALCULATES  THE  DEFIIITE  IITEGRAL  OF  ** 
C  COSIIE(X)  RAISED  TO  AI  IITEGRAL  POWER.  FORMULAS  « 
C  ARE  TAKER  FROM  GRADSHTEYI  ft  RYZHIK.  ** 
C  ** 


O ******************************************* ************************** 

C  »»  CALLED  BY:  IMIABC 

C  »»  SUBROUTIIES  CALLED:  IOIE 

C  »»  FUHCTIOHS  CALLED:  BIIOMIAL 

C 

C  »»  IITERHAL  VARIABLES  «« 

C  I  =  DO-LOOP  IIDEX 

C 

C  »»  DATA  IIPUT  FROM  CALLIIG  ROUTIIE  «« 

C  I  =  POWER  TO  WHICH  COSIIE  IS  RAISED 
C  A  =  LOWER  LIMIT  OF  IITEGRATIOH 

C  B  =  UPPER  LIMIT  OF  IITEGRATIOH 

C 

O********************************************************************* 

REAL*8  FUICTIOI  IITEGRAL_COSIIE(I,  A,  B) 

C 

IHTEGER  I,  I 
REAL*8  A,  B,  BIIOMIAL 
C 

C  IF  1=0,  IITEGRAL  IS  TRIVIALLY  EASY  TO  EVALUATE 
C 

IF  (I.Eq.O)  THE! 

IITEGRAL_COSIIE=B-A 

C 

C  IF  I  IS  EVER,  USE  GRADSHTEYI  ft  RYZHIK’ S  EQUATIOI  2.513.3 
C 

ELSEIF  (M0D(I,2).IE.l)  THEI 

IITEGRAL_COSIIE= (BIIOMIAL (I ,  I/2)/(2 . **I))*(B-A) 

DO  10  1=0,  (I/2)-l 

10  IITEGRAL_C0SIIE=IITEGRAL_C0SIIE+(0.S**(I-1))* 

ft  (BII0MIAL(I,I)/(I-2. *1))* 

ft  (DSII((H-2.*I)*B)-DSII((I-2.*I)*A)) 

C 

C  IF  I  IS  ODD,  USE  GRADSHTEYI  ft  RYZHIK ’S  EQUATIOI  2.513.4 
C 

ELSE 

IITEGRAL_C0SIIE=0 . DO 
DO  20  1=0,  (I-l)/2 

20  IITEGRAL_C0SIIE=IITEGRAL_C0SIIE+(0. 5**(I-1) )* 

ft  (BIIOMIAL (I, I)/ (1-2. *1))* 
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(DSII((I-2.*I)*B)-DSII((I-2.*I)*A)) 


* 

ENDIF 

RETURI 

END 
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5.7  Function  SINC 


C********************************************************************** 


c  ** 

C  »»  FUHCTIOH  SIHC(X)  ««  ** 
C  ** 
C  THIS  FU1CTIOI  RETURIS  THE  SIIC  FUICTIOI,  OR  SII(X)/X,  OF  ** 
C  A  REAL  HUMBER,  X.  ** 
C  ** 


Q* **********************************************************  *********** 

REAL*8  FUHCTIOH  SIHC(X) 

C 

REAL*8  X 

C 

C  CALCULATE  SIHC,  USIIG  THE  LIMITIHG  FORM  IF  X  IS  LESS  THAH  O.OOOl 
C 

SIHC=1 .DO 

IF  (ABS(X) .GT. 0.0001)  SIHC=DSII(X)/X 
C 

C  ELIMIHATE  ROUID-OFF  ERROR  IF  HEAR  A  ZERO  CROSSIIG 
C 

IF  ((ABS(SIHC).LT.1.D-12). AHD.CX.LT. 100))  SIHC=O.DO 
C 

RETURH 

EVD 
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3.8  Function  CONJUGATE 


C********************************************************************** 
C  ** 

C  »»  FUNCTION  COS JUGATE(X)  ««  ** 

C  ** 

C  THIS  FUSCTIO!  RETURNS  THE  COMPLEX  COIJUGATE  OF  A  COMPLEX  ** 

C  NUMBER,  X.  ** 

C  ** 

C********************************************************************** 
COMPLEX* 16  FUHCTIOH  COHJUGATE(X) 

C 

C0MPLEX*16  X 
C 

C  CALCULATE  CONJUGATE 
C 

COM  JUG ATE=REAL ( X )  +  (0 . ,-1 . )*DIMAG(X) 

C 

RETURN 

END 
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S.9  Function  ARG 


C*** ******************************************************************* 


C  ** 
C  »»  FUHCTIOH  ARG(X)  ««  ** 
C  ** 
C  THIS  FUICTIOI  RETURHS  THE  PHASE,  OR  ARGUMEHT,  OF  A  COMPLEX  ** 
C  HUMBER,  X.  ** 
C  ** 


C********************************************************************** 
REAL*8  FUHCTIOH  ARG(X) 

C 

COMPLEX* 16  X 
C 

C  IF  |X|=0  THE H  RETURH  ZERO  AID  FLAG  THE  ERROR 
C 

ARG=O.DO 

IF  (X.EQ. (0. ,0.))  THEH 

WRITE  (6,*)  'Undefined  phase  in  subroutine  ARG' 

RETURH 

EHDIF 

C 

C  CALCULATE  ARG(X) 

C 

ARG=DATAH2(DIMAG(X) , DREAL(X) ) 

C 

RETURH 

EHD 
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